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ABSTRACT 

Techniques for the identification of limit cycles in nonlinear 
control systems rely heavily upon .describing function analysis. 

Most methods are restricted to single nonlinearity systems' and 
consider only the fundamental component of the Fourier series 
representation of the nonlinearity- output. No previously available 
technique considers multiple harmonic components in systems contain- 
ing several nonlinearities and multiple forward paths. 

This report presents a method of analysis which identifies 
limit cycles in autonomous systems with multiple nonlinearities in 
multiple forward paths. The nonlinearities may be any time- and 
frequency- invariant single input, single output type which can be 
adequately repr esented in piecewise-linear form. In addition, 

several harmonic components may be retained in the procedure. 

/ 

Included in the report is the FORTRAN code for implementing 
the Harmonic Balance Algorithm (HBA) developed in the research. 

Using this code, limit cycles are identified in multi-path, 
multi-nonlinearity systems while retaining the effects of several 


harmonic components. 
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CHAPTER I 
INTRODUCTION 


The science of control system analysis has experienced rapid 
development in the years since WWII. Many techniques for 

* 

determining the stability of linear feedback control systems have 
evolved, some of which give a simple indication where others give 
additional information such as the margin of stability. These 
methods are described in detail in many texts, such as Truxal [1], 
Dorf [2], and Kuo [3], and are summarized here. 

In a system such as that shown in Figure 1.1, a transfer 
function is the ratio of the output and input of a block. The 
transfer function of the upper block is thus G(s) = C(s)/E(s) 
while that of lower block is H(s) = V(s)/C(s) . By algebraic 
manipulation the input-output relationship of the system is 
obtained as 

C(s) _ G(s) ,, .. 

R(s) 1 + G(s)H(s) * K } 

The denominator of the above equation when set equal to zero is 
referred to as the characteristic equation; it determines the 
mathematical characteristic of the transient response of the system; 


* A linear system is one for which the superposition principle 
holds. Superposition implies the properties of additivity and 
homogeneity with respect to both input and initial conditions. 
Additivity is illustrated in the following manner: if input rj 
produces output yi , and input r 2 produces output yi , then 
input ri + r 2 will produce output yi + y2 • Homogeneity is 
illustrated thusly: if input ri produces output y 1 , then 

input kri will produce output kyj . A linear system is stable 
if its output tends to a finite value for any finite input. 
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A relatively simple mechanistic procedure for determining sta- 
bility from the unfactored polynomial in s , which forms the charac- 
teristic equation, is the Routh test. This test, which originated 
about 1880 , can be used to determine the number of the roots with 
positive real parts. The existence of roots with positive real parts 
indicates an unstable system.. 

In addition several graphical techniques are widely used in . 
testing for stability. These, in general, provide information re- 

i * -;’ u 

garding the relative stability of the system! Included are the 
Nyquist Criterion, the Bode Diagram, and the Root Locus Plot. 

* i * 

There are cases where the transfer function G (or H) cannot 

be readily determined as a polynomial in s . When the polynomial 

* * 

is not available, the Routh test cannot be applied. The Nyquist 
stability criterion . can be used in these cases if the frequency 
response characteristics of G(s)H(s) can be determined. 

Note that for a linear system to be stable, all zeroes of 
1 + G(s)H(s) must lie in the left half of the s-plane. The Nyquist 
criterion determines if this is the case by examining a plot of 
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G(s)H(s) as s is allowed to vary along a closed contour in the 
s-plane. This contour may be, for example, the imaginary axis of 
the s-plane together with a large semicircle enclosing the right 
half plane. An indication of stability can be obtained by examining 
the number of resulting .clockwise encirclements of the point (-1,0) 
in the G(jto)H(jto) plane. The number of encirclements is given by 

N = Z - P (1-2) 

=* the algebraic sum of encirclements of the (-1,0) 
point where clockwise encirclements are assumed 
to be positive. 

= the number of zeroes of 1 + GH in the right 
half of the s-plane 

= the number of poles of GH in the right half 
of the s-plane 

If there are any net clockwise encirclements of the (-1,0) 
point, there are zeroes of 1 + GH in the right half of the 
s-plane, and the system is unstable. The number of right half 
plane zeroes is given by the sum of N and P . In many cases the 
open loop system is known to be stable, i.e, there are no poles of 
GH in the right half s-plane. (Note that if the frequency response 
characteristics of GH are determined by measurement, GH must be 
stable.) Thus N must equal zero for the system to be stable. 

Nyquist plots are normally made using polar notation. The 
magnitude and phase of GH(s) is calculated and plotted for 

4 

various values of to when s is replaced by jto . Thus frequency 


where N 


Z 


P 
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information is not directly available from the plot; it must be 
obtained by estimation for points on the plot which, are located 
between calculated points. 

The Bode method, developed by Hendrik Bode of the Bell 
Telephone Laboratories in the 1930's, is a graphical technique 
which provides frequency information directly. This technique 
utilizes two graphs: the magnitude of GH(joj) and the phase of 

GH(jw) , each plotted as a function of m . For ease in construc- 
tion and interpretation, logarithmic scales are commonly used for 
both the frequency and gain axis. For well-behaved systems an 
indication of stability is obtained by an examination of the gain 

at the frequency which results in 180° of -phase, or the phase at 

* 

the frequency which produces a gain of unity. In the former case 
the gain should be less than unity while in the latter case the . 
phase lag should be less than 180°. 

The final method of linear system analysis to be discussed 
was developed by Evans [4] in 1948. This technique, referred to as 
the Root-Locus method, is a graphical representation of the s-plane 
loci of the variation of the poles of the closed-loop system func- 
tion with changes in open-loop gain. [5] The poles and zeroes of 
the open loop system are located in the s-plane; then the loci of 
the closed-loop poles are plotted as some parameter is varied 
(usually open loop gain) . If these loci are totally contained 

* Well-behaved systems are defined as those which do not have 
any open loop poles in the right half plane and only have one 180° 
crossing and/or one unity gain crossing. For systems which are not 
well-behaved, stability is best studied with the Nyquist criterion. 



5 


within the left half of the s-plane, the system is stable. Should 
loci cross into the right half plane, the system will be unstable 
for the values of the parameter that result in that part of the 
locus in the right half plane. 

Unfortunately the well developed techniques of linear systems 
analysis do not in general extend to nonlinear' systems. Some 
authorities [6] hold that a general method of synthesis is impossible 
since a nonlinear system is one whose response -is dependent upon 
its input. Nevertheless, techniques have been developed which apply 
to certain types or classes of nonlinear systems or which give 

limited information about system performance. Several of these 

« 

techniques are discussed briefly here and are presented in detail 
in texts such as Hsu and Meyer [7] and those previously mentioned. 

There' are several unusual characteristics of .nonlinear systems 
which are of interest. These are defined here and will be referred 
to throughout this_ discussion. 

1. LIMIT CYCLE: Limit cycles are oscillations of fixed 

amplitude and period which can occur in nonlinear systems. 
There are two types of limit cycles, stable and unstable. 
When a system is operating in a stable limit cycle, any 
small disturbance results eventually in a convergence to 
the original limit cycle. On the other hand, any small 
disturbance to a system operating in an unstable limit 
cycle results in eventual operation at .some other point. 


A system may exhibit behavioral modes that include both 
stable and unstable limit cycles. 
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2. HYSTERESIS: Hysteresis is a nonlinear phenomenon wherein 

the output depends not only upon the value of the input, 
but also upon its history. A familiar example is the 
conventional magnetization curve. 

Many nonlinear systems, in which the deviation from linearity 
is not too great, may be examined' using linear techniques. In these 
quasi-linear systems, the nonlinearity is replaced by a linear 
approximation of its characteristics. While the resulting analysis 
is approximate, this approach is useful under the concept that an 
approximate answer is better than no answer at all. Although 
practical systems are almost never purely linear they are often 
quasi-linear; this is precisely why techniques for analyzing linear 
systems produce such good results. 

The describing function is an attempt to extend transfer 
function analysis techniques developed for linear systems to the 
analysis of quasi-linear systems. The normal form of this technique 
is a frequency-response method which is based on an analysis which 
neglects the effects of all harmonics except the . fundamental . It 
is therefore most successful when applied to systems which contain 
components that are low-pass. 

A single-loop system with one nonlinearity is shown in 
Figure 1.2. Since the describing function technique applies 
directly only to autonomous systems, the system is more properly 


* An extension of the normal describing function method is the 
dual input describing function, which retains two harmonics or the 
fundamental and dc term. The dual input describing function is not 
considered here. 
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Figure 1.2. 


Single-loop system with nonlinearity N 



Figure 1.3. Autonomous system with nonlinearity N 
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considered in the forms of Figure 1.3. It is assumed that the 
linear element transfer function G is low-pass and allows only 
the fundamental component of m to pass. Therefore c , and thus 
e are sinusoidal and may be written as 

e = E sin cot (1-3) 

The output m of the nonlinearity may be expressed as 

m = K(E) + f ,(e) (1-4) 

a 

1 

where K is a- quasi-linear equivalent gain of the nonlinearity and 
f^ is a distortion term. 

The distortion term f.^ .may be minimized in ‘a mean-square 
sense by choosing K as the normalized coefficient of the funda- 
mental of the output wave. Thus with an input e as given in 
equation (1-3) the Fourier series representation of the output will be 

m <= hj sin cot + cos cot + ... (1-5) 


The describing function is defined, in cartesian- coordinates, as 


K 




(1-6) 


In general, for symmetric single-valued nonlinearities there will 
be no phase shift in output fundamental; the imaginary part will 
be zero. There will be phase associated with double- valued, or 
hysteresis, functions and with nonsymmetric nonlinearities. 
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Applying the describing function approximation to N yields 

m - -Kc (1-7) 

If a limit cycle exists only the fundamental component is retained 
in this approximation. The steady-state form of the linear transfer 
function can thus be used; that is 

G( « ' a ~ 8) 

If indeed a limit cycle does exist, (1-7) and (1-8) must be 
satisfied simultaneously., giving 

G(j «) - ~~ (1-9) 

This relationship is the basis of a graphical analysis technique 
making use of the Nyquist criterion ♦ 

First, on a polar graph a Nyquist, plot of the linear system is 

+ 4 » ' 

made as illustrated by G(jui) . in Figure 1-4. Next, a plot of - — . 

K 

is made as a function’ of E on the same polar graph. Oscillations 

I • 

' -I 

are indicated by intersections of the G(jaj) and the ' - — curves 

K 

as indicated by (1-9) . 

If the describing function contains no phase shift, i. e. , the 
imaginary part is zero, the plot of — ^ lies along the negative 
real axis, as illustrated by curve kx in Figure 1.4. If the 
nonlinearity can cause phase shift between its input and output a 
curve as illustrated by k 2 can be expected. If the G(jto) curve 
and the k curve intersect, limit cycles will exist with approxi- 
mate frequencies being those at the intersections. The values of 
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Figure 1.4. Nyquist-describing function plot of 
system of figure 1.3 
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E at which the intersections occur are the approximate amplitudes 
of the limit cycles. 

Whereas the describing function technique is an adaptation of 
linear system analysis methods to the study of nonlinear systems, 
the phase plane method involves a radically different approach. 

This analysis is concerned with the characteristics of the solution 
of the differential equation 

x + a(x,x)x + b(x,x)x = 0 (1-10) 

The terms a(x,x) and b(x,x) are functions of the signal x and 
its derivative x . The derivative x is plotted as a function of 
x producing a phase plane portrait. The initial conditions x(0) 
and x(0) , for time t = 0 , define a point in the phase plane; 
the behavior of the system at all later times is identified by the 
curve whose path passes through this defined point. If as time 
increases the path tends to infinity, the system is unstable; if 
the path tends “toward the origin, the system comes to rest and is 
stable. There is also the possibility that the path may close on 
itself; this indicates the existence of a limit cycle. 

Around the turn of the century, Russian mathematician A. M. 
Lyapunov presented a useful stability criterion which has found 
extensive application in the past 20 years. This approach examines 
the response of a system at an equilibrium point when it is dis- 
turbed slightly. An equilibrium point x g is stable in the sense 


* An equilibrium point is a state in which the system is in 
equilibrium with no tendency to move unless driven. 
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of Lyapunov if, for a given e -> 0 , there exists a 6 > 0 such 
that when the initial disturbing condition x ' satisfies 

II* 0 -* e II <s (i-ii) 

then 

||x(t)-x |j' <e for all' t >_ 0 . (1-12) 

A more restrictive sense of stability requires that the system 
return to the equilibrium point after disturbance. This condition 
(termed asymptotic stability*) is met if, in addition to the previous 
requirements for stability (Equations 1-11 and 1-12) , 

lim x(t) = x (1-13) 

t-*» e 

for I |x q - xj I < 6 

These conditions can be shown to be true if a function V(x) , 
called a Lyapunov function, can be found which is positive every- 
where (except possibly zero at x^) , and whose derivative is negative 
everywhere except at the equilibrium point. The function V(x) is 
analogous to (but in general not equal to) the energy stored in 
the system. If it is everywhere positive but decreasing, the 
system is continually settling back toward its equilibrium point 
and is asymptotically stable. Note that any function V(x) which 
satisfies these conditions proves asympototic stability. It should 
be noted that V(x) is not necessarily unique. 

The techniques presented here represent a cross-section of 
methods useful in control system analysis; it is far from an exhaus- 
tive listing. Several other methods are discussed in the previously 
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referenced texts and in the current literature. Generally these 
other methods either are restricted to a select class of systems, or 
contribute no information about stability which could not be obtained 
using the described methods. . However, since each method has certain 
limitations and strengths, one finds it useful to apply several tech- 
niques when carefully examining a particular system. Insight to 
system characteristics is frequently gained from one analysis which 
is not evident using other techniques. 



CHAPTER II 


THE PROBLEM 

The analysis of nonlinear feedback control systems is generally 
a complex task requiring a variety of tools. Sometimes a complete 
characterization of the system is not necessary; identification of 
limit cycles often yields sufficient information. ‘ The methods 
discussed previously are quite useful and valuable; however, they 
each suffer from certain limitations. 

The technique of Lyapunov offers an important advantage in 
that stability of .systems can be determined without having to solve 
the system equations. The method does require the identification 
of a Lyapunov function, however, which merely provides sufficient 
conditions for stability. If a function of the required type cannot 
be found for a given system it does not mean that the system is 
unstable. It means simply that the attempts to establish stability 
have failed. If, on the other hand, stability is established, the 
resulting region of stability may be more restrictive than necessary; 
stability may exist in a region larger than that specified by the 
method . 

In addition, no systematic procedure exists for evaluating an 
arbitrary function of. several variables. As a result, test funotions 
are frequently restricted to the quadratic form, greatly limiting 
one’s choice of test .functions. 
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The graphical phase-plane method, in addition to providing 
stability information,' also provides some insight into system 
characteristics. Limit cycles and stable equilibrium points are 
readily identified and some indication of system response is 
obtained. The technique is unfortunately limited to second-order 
systems since higher-order derivatives cannot be conveniently 
displayed on a two dimensional plot, (Knowledge of the variable 
and its first derivative would not completely define a higher-order 
system.) Also, as with all graphical techniques, the accuracy of 
the results is at least partially dependent upon the care used in 
the graphical construction methods. These methods are often tedious 
and somewhat prone to error. 

Perhaps one of the most widely used techniques for determining 
limit cycles in practical control systems, particularly those of 
higher-order; is the classical describing function method. This 
method is based on an analysis which assumes that the signals in 

the system can be completely .characterized by a truncated Fourier 

1 

series which neglects the effect of most harmonics. With only a 
single nonlinearity, this assumption may yield incorrect or 
misleading conclusions. [8] In considering multiple nonlinearities, 
questionable results are even more likely. 

The graphical approach normally associated with the classical 
describing function is not suitable when considering the effects of 
retaining more than the fundamental component of the truncated 
Fourier series. An extension of the describing function technique 
allows the consideration of two components. This is accomplished by 
use of the dual-input describing function, but at the expense of 
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increased complexity of the analysis technique. In addition, these 
techniques are limited to the consideration of single path, single 
nonlinearity system^, .although’ (as- will be -seen later) this restric- 
tion can sometimes be relaxed. 

Gelb [9] developed a method of analysis of systems of multiple 
identical nonlinear-linear combinations. He observed that in a 
system such as that of Figure 2.1, for nonlinear operation on a limit 
cycle to exist, x^ = except for a time phase difference. The 
method is applicable only if the linear-nonlinear combinations are 
identical. „ 

, Gronner [10] developed a method of analysis for systems having 
cascaded nonlinearities. The cascaded nonlinearities are first 
reduced to a composite; the describing function application of the 
Nyquist plot is used to identify any limit cycles. The limitations 
usually associated with the describing function technique apply to 
the resulting composite. 

Gran and Rimer [11] developed a modified Root-Locus technique 
using the describing function representation of any nonlinearity to 
give closed-loop analysis. This technique may be used to determine 
the existence of any limit cycles and to obtain stability information. 
The method works for both amplitude and frequency sensitive non- 
linearities, and for multiple nonlinearities separated by dynamics. 
Since system dynamics are displayed in pole-zero configuration, the 
technique is useful for system synthesis. Because the describing 
function is used, the previously enumerated limitations apply. 
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Jud [12] devised a method using the describing function 
technique to determine limit cycle conditions for any number of 
parallel nonlinear elements as long as signals reaching each 
nonlinear element have passed through a low pass filter. The 
implication is that each nonlinearity is preceeded by a low pass 
linear element. Thus, only a fundamental component of frequency 
reaches the nonlinearity, and the describing function produces 
only a fundamental component on the output of the nonlinearity. 

Negoescu and Sebastian [13] introduced a method for determining 
self-oscillation in a nonlinear system which contains two nonlinear- 
ities separated by a linear device. The nonlinearities appear in 
parallel feedforward paths so constructed that in all paths each 
nonlinear element is separate from the other nonlinear element by 
a linear element, as shown in Figure 2.2. System equations .are 
written, then examined for solutions. The existence of roots with 
positive real-parts indicates oscillations. 

The most general method of analysis applicable to systems 

with linear and nonlinear elements was developed by Davison and 

Constantinescu [14]. While the technique is based on a single-loop 

system as shown in Figure 2.3, the authors have extended its use 

to include parallel systems of the type depicted in Figure 2.4. 

The normal describing function approximations are assumed, i.e., 

* 

the linear elements are sufficiently low pass that, if a limit 
cycle exists, only the fundamental frequency component need be 


considered. 




Figure 2.2. Parallel multiple nonlinearity system 
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Under these conditions, the inputs' to the nonlinearities take 
the form of 

ej(t) = Ai sin tot 

e2(t) = A2 sin (tut + ©2) ( 2 - 1 ) 

♦ 

» ' 

e (t) = A sin (cot +. 0 ) 

n n 

Since the describing functiqn may be considered a linear gain 

1 

whose magnitude M and phase <J) are dependent upon the input mag- 
nitude, the nonlinearity outputs may be written as 

m‘l(t) = MiAisin({ot+<f>i) (2-2) 

m 2 (t) = M 2 A 2 s in (wt+0 2 + 4>2 ) 

« 

m (t) = M A sin(wt+0 +<f> ) 

n' n n n n 

If the magnitude of the sinusoidal input to the nonlinearity is 
fixed, conventional linear analysis techniques may be applied to 
the system. 

For the limit cycles to exist 

S’ G (j»*) G d (A ) = -1+jO (2-3) 

i=l 

where w* is the frequency of oscillation of the limit cycle, 
Gdi(A^) is the magnitude and phase of the linear gain equivalent of 
the ith nonlinear element and is the gain of the i fc ^ linear 


element . 
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In addition, each linear element gain may be expressed in 

4 

terms of the nonlinear gains and inputs. .Thus 


G^j®*) 

G n (j®*) 


fe i+l 

e i G d i (A i' ) 



i ^ 1 , 2 , • . * , n -1 


(2-4) 


When combined with . (2-1) , (2-4) may be written 

A 2 = Aj | Gi (j (o*) I | G dl (A x ) 

A 3 = A 2 !G 2 (jco A )||G d 2 (A 2 ) (2_5) 

A limit cycle exists if there is. an Aj, A 2 , A 3 , ... and w which 
simultaneously satisfies both ('2-3)and £2-5) . A graphical method may be 
used to determine any such solutions. The following procedure is 
used. 


1 . 


2 . 


3. 


Plot jj G^(jo)) as a function of w . 
i=l 

Choose an arbitrary value for Ai 
n n 

Plot H • G.(jw) n G, (Aj_) as in Figure 2v?5, 
i=l 1 i=l 1 


using equa- 


tion 2 r *5 and step ( 1 ) . 

4. Determine possible limit cycle frequencies by noting the 
frequencies at which the plot of step (3) intersects the 
negative real axis. 


**If the number of memory type nonlinearities does not exceed 
one, a somewhat simplified procedure may be used: See Davison, et. 

al. f 141 for details. 
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5*. Determine the magnitude k of the plot of step (3) at 
frequencies determined in step (4) . 

6. Repeat steps (2)-(5) for different values of . Plot 
the graph of k versus A as shown in Figure 2*6* The 
intersection with the k = 1 line identifies limit cycle 
frequency-amplitude values* 

7. Using Equation 2-5, determine the remaining A^ values* 

The system of Figure 2.4 (note that there is only one nonlin- 
earity in each parallel path) may be analyzed in a similar fashion 
with the following modifications in notation: 

The input e is given by 

e = A s^n a) t (2-6) 

and the nonlinear element inputs are 

e^.(t) = A^ sin( o)t-H>^) (2-7) 

where 

A i = | H ± (j to) | A (2-8) 

4>i = L (HjOoO i = l,2,...n (2-9) 

The conditions for the existence of limit cycles are Equation 
2-8 and 

l H.(joj) G (j») G d (A ) = -1 + j0 . (2-10) 

i=l 1 1 i 1 

The graphical procedure is as follows: 


1) Choose A arbitrarily. 
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n 

2) Plot £ H. (joj) G.Cjw). G<j (A.) as' a function of w , 

m n. n 1 

1=1 1 


3) 


4) 


making use of Equation 2-8.; 

Determine possible limit cycle frequencies from the inter- 
section of the plot of step (2) and the negative real axis. 

i n I 


Determine k = 


I H. (jio) G (jco) G (jo)) G(j (A i ) 
i=l 1 


at the frequencies determined in step (3) . 

5) Repeat (l)-(4) for different values of A . Plot k 
versus A, noting the points of intersection with the 
k = 1 line. 

6) Use (2-8) to determine corresponding frequencies of 
oscillation and input amplitude A^ for the limit cycles. 

While the Davison-Constantinescu method provides a powerful 
technique for the analysis of systems with multiple nonlinearities, 
it suffers from the normal restrictions of all techniques using the 
describing function approximation. This approximation is based on 
the assumption that elements of the system are assumed to be low- 
pass. Thus the input to each nonlinearity is accurately represented 
by a sinusoid. If harmonic components are present at the nonlin- 
earity inputs, results obtained by any method based on describing 
function analysis are suspect. 

Thus there is a need for a method of analysis which relaxes 
the describing function restriction. The work presented in the 
following chapters removes the low-pass requirement of the system. 
The nonlinearity inputs may therefore contain a dc component 



28 


plus several harmonics, as well. as the fundamental frequency. In 
addition, parallel paths may contain multiple nonlinearities. 



CHAPTER III 


THE SOLUTION 

In the last chapter the difficulties and problems associated 
with using describing function techniques to determine harmonic bal- 
ance in feedback control systems were delineated. ■ In this ^chapter an 
iterative approach for .analyzing a general', class.- of - nonlinear -feed-- 
back control systems for harmonic balance is developed, along with 
the necessary mathematical tools for implementation on the digital 
computer . 

3.1 Harmonic Balance 

A single loop time-invariant nonlinear feedback control system 
is shown in Figure 3.1. Under conditions of harmonic balance the 
signal at each point in the loop will be periodic- and of the same 
period. Thus, each signal, can be represented as a Fourier series, 
i.e. , 

CO 

e (t) = a + 7 a cos kco t + b, sin kto t (3-1) 

° k=l k ok 

sin km t 
o 


CO 

c (t) = A . + J A^ 'COS kw t + B, sin km t 
0 k=l ^ ok o 


u (t) = a + l a, cos kiivt + $ 
° k-i k 0 
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Figure 3*. 1 ♦ Single loop nonlinear, 
time-invariant control system. 

The necessary and sufficient condition for harmonic balance is 

e(t) - -c(t) . (3-2) 

Therefore, each term in the series for e(t) must be equal to the 
corresponding term for -c(t) . This requires that 

a = -A (3-3). 

o o 

a^ = k — 1,2,3,... 



In processing the signal the linear element may alter the 
magnitude and phase of the various harmonics as a function of their 
corresponding frequencies, but no new harmonics will, be generated. 
The nonlinear element, however, can significantly alter the 
harmonic content of the signal. The establishment of harmonic 
balance depends upon the selection of the magnitudes and frequent 
cies of the signal such that the coefficient equalities in C3~3)‘ are 
satisfied. This requires the solution of an infinite number of 
simultaneous, nonlinear, algebraic equations in an infinite number 
of unknoxms (obviously an ■itnpcu>6 &') . 





31 


Fortunately the problem can be somewhat simplified. In most 
practical systems the linear or nonlinear element can be assumed to 
be band limited, which results' in a finite number of terms in the 
Fourier series representation of the output signals.* 

t 

The problem then becomes one involving only a finite number of 
equations in a finite number of unknowns, which can be solved using 
iterative procedures. Such procedures examine the output produced 
from an initial guess of the input, then refine the- input to-, 
produce an output that more nearly satisfies (3-3) . This process 
is repeated until the output matches the input within a desired 
tolerance. 

3. 2 General System Considered 

The system of Figure 3.1 can be expanded to include several- 
paths with multiple nonlinearities in each path, as shown in 
Figure 3.2. In order for harmonic balance to exist 

l - c = - e (3-4) 

£=1 

for all components of c and e. That is, the dc component of e 
must he equal to the negative sum of the dc components of the Y f s; 
a similar condition is true for each component of the Fourier series 
representation of the Y*s and e. 


* It should be noted that the describing function technique 
solves the harmonic balance problem in this way. In the describing 
function case only one or two terms are retained in the .Fourier 
series. 
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The goal of the work presented here is to develop a computer 

■v 

code to determine harmonic balance in nonlinear feedback control 
systems with the general configuration shown in Figure 3.2. The 
technique that is developed on the ensuing pages is designed to 
contend with several Fourier series terms, in the establishment of A 
harmonic balance (thus making it superior to the -describing 
function methods mentioned in Chapter 2) . 

In the development of any technique of this type some limiting 
assumption are necessary. In this case they are: 

1) The system can be separated into linear and nonlinear com- 
ponents. 

2) No input signal is applied. 

3) The nonlinear element is ’time- and frequency- invariant. 

'4) The input /output characteristics of the nonlinear element 

can be approximated by piecewise-linear segments. 

5) The nonlinear elements are single input/single output. 

6) The linear elements are relatively low-pass. 

3.3 The Harmonic Balance Algorithm (HBA) 

The flow chart in Figure 3.3 shows the major steps which must 
he undertaken to solve iteratively the harmonic balance problem. 
Using an initial guess of a limit cycle (including frequency and 


*It will be assumed that no frequency component higher than the 
highest harmonic retained in the analysis will be passed by an eln- 
ment in the system. 
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Figure 3.3. Flow chart for HBA 








35 


amplitude components) the system output is determined. For the 
present a "No" and "Yes" are, respectively, assumed for the two 
decision blocks. An input change vector is then determined so as 
to reduce the difference between the input and the output components. 
Using this change vector, a step size is selected, and the input 
values are incremented. 

This new input is utilized in computing a new output. If the 
new output represents an improvement, the process continues through 
another cycle. If not, the selected step r size, is reduced and new 
input values are determined. This input is then used to recalculate 
the output, initiating a new cycle, 

. Each time new output values are computed they are compared to 
the input values to determine if another iteration is needed. The 
process continues until output an<} input agree within the desired 
tolerance, at which time the procedure terminates. 

The details of the process together with the mathematical basis 
for the procedure are presented in the following sections. The 
FORTRAN code that implements the algorithm is presented in Appendix B. 
Directions for using the HBA and q sample output are presented, 
respectively, in the Appendices A and C. 

3.4' Determination of System Output 

Harmonic balance is established by the necessary and sufficient 


condition that 
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(3-5) 

where the subscripts correspond to the Fourier components of the 
variable. The procedure described in the last section for accom- 
plishing', this equality involves several distinct operations' (the 
efficient performance of which ate central to the development of 
a pragmatic analysis tool) . 

First, the computation of the Fourier series of the output of 
nonlinearities when the input is a truncated Fourier series is 
necessary. At the onset of this research effort several methods 
were' investigated: a Simpson’s numerical integration technique, a 

Kalman filter technique, and a method devised by Braswell [15] in 
his Relaxation-Based Algorithm. As a result of testing, Braswell's 
technique was ultimately chosen on the basis of comparable accuracy 
with somewhat faster execution time. 

The Braswell procedure uses the secant method (also known as 
Regula falsi, or method of false postion) to determine the times 

{t.} which identify the nonlinear breakpoint crossings by the 

1 i=l 

input signal. These times are illustrated in Figure 3.4. Each 
breakpoint is identified by a coordinate pair consisting of the 
value of the input and the resulting output value at the junction 
of two of the linear-piecewise segments. The slope of each segment 



is also specified. 
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Once-:'each breakpoint time -is identified the proper breakpoint 
coordinate and slope are linked to it. With this information, then 


t. 

I f 1 


d m T B i • 

x 1 t.y-i 


e “j ffi Vdt 


I H 

+ T .^i l C n 
1=1 . T 
rt =-N 


l 

f 


e j (n ‘ m) Vdt 


V 1 


(3-6) 


is used to calculate the resulting complex Fourier series coeffi- 
cients of the nonlinearity output. In (3-6), d^ is the Fourier 
coefficient of the harmonic term and is the slope of the 

piecewise- linear segment used on the interval [t^, t^j^] . Also 

B. 83 u. - -K.e. - ; this relationship identifies the ordinate 
i l-l i l-l r 

intercept of the piecewise- linear segment. 

In processing the signal, linear elements may alter the 
magnitude and phase of the various harmonics, but no new components 
will be created. Since the principle of superposition applies, the 
linear element output can be obtained by multiplying each harmonic 
component of the nonlinearity oufput by the gain and phase of the 
linear element at the harmonic frequency. Thus, 

Y = G(jmw ) d m — 0,1,2, . . . ,'M (3-7) 

where the m represents the zeroth, first, , and Hth harmonics. 

The technique is applied to each path of Figure 3.2 and the 
results are summed to produce the system output. That is, the out- 
put of the first linear element in path one is taken as the input to 
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nonlinear element two to obtain the output of linear element two. 
This, in turn, is used at the input of nonlinearity three to deter- 
mine'fhe output of linearity three. This process is continued until 
all linear-nonlinear pairs of path one are considered. ' Then the 
original input is assumed at nonlinear element one of path two and 
the path two output is determined. The process is repeated until 
all path outputs are found and their sum is obtained. 

This procedure allows one tq determine the output which results 
for a given input. Some means are .necessary to determine the changes 
in input which would result in the solution of the relationship 

c - e = 0 , (3-8) 

n n 

giving harmonic balance. This becomes a zero finding problem in 
which the variables are frequency, wi, and the Fourier coefficients 
of the input (and the output) signals. 

In actuality each input consists of a Fourier series with an 
infinite number of terms. Thus, it is necessary to solve an infi- 
nite set of algebraic equations in an infinite number of unknowns. 
However, since the linear elements are assumed to be relatively low 
pass, only a finite number of terms will be significant. The 
resulting nonlinear simultaneous equations are of such complexity 

that they do not lend themselves to an exact analytical solution. 

A numerical solution technique base on the Newton-Raphson method is 


devised in the next section. 
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3.5 Determination of Trial Solutions 

Another operation in Figure 3.3 is the computation of the in- 
put directional change vector. The Newton-Rap hson method is an 
iterative technique that can be used for solving a system of simul- 
taneous nonlinear equations [16] (such as those found in this pro- 
blem) . The method automatically selects the increment to be used 
in determining the next iteration input and as a rule exhibits a 
rapid rate of convergence (if it converges) . The Newton-Raphson 

method is selected for the solution tool in this work. However, it 

* 

is modified to improve convergence. 

The technique will first be illustrated for a system of two 
equations, i.e. 


f 1 (x 1 , x 2 ) = 0 (3-9) 

f 2 (Xi » x 2^ = 0 


An initial guess is made for the unknowns x^ and x 2 . 

Small increments Ax^ and Ax^ are assumed to give the correct solu- 
tion to (3-9). That is. 


f l ^ X 1 + Ax l ’ X 2 + Ax 2^ = ° 
f 2 (x° + Ax^ , x° + Ax 2 ) = 0 


(3-10) 


where x° , x° , Ax^ , and x 2 are respectively the initial guesses 
and increments. 


*This results in essence in a gradient procedure. 
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Expanding Equation 3-10 in a Taylor series around the initial 
guess yields 


(3-11) 


^1^ X 1 * X 2^ + ^ x l ^ iL ^ ( 1 ) + ... = 0 

8x l 2 8x 2 


n rt 3f 9 3f, 

f 2 (x 1 , X 2 ) + ^ ^ > +& 2 C.a^ 


) + . . . =0 


where the partial derivatives ape evaluated at the initial guesses, 
x^ and . 

If the increments are assumed to be small., higher-order terms 
of Equation (3-11) are negligible. Equation (3-11) may thus be 
written in matrix-vector form as 


where 


A 

J = 


• j L'x; 

* 0 

3f i 

3f l 

3x i 

Sx 2 

3f 2 

3f 2 

3X 1 

3x 2 


(3-12) 


(3-13) 


(the Jacobian matrix of the functions f^ and f 2 with respect to 


x^ and x 2 ) , 


F = 


f • x 2 > 


f (x° x° ) 
X 2 K 1 * 2 } 


(3-14) 



Solving Equation (3-12) for AX yields 


AX * -[J] 1 F 


(3-16) 


The components of AX are Ax^ pnd Ax£ . When these values are 
added to the initial estimates x° and x° an approximate solution 
of (3-10) is determined. However, since they -are approximate, the 
solution is not exact and additional refinement in the values of 
the x's are needed. This is obtained by using the values of 
x° + Ax^ and x° + AX 2 as new initial guesses of x^ and X 2 
and performing the operation again. The process is repeated, 
obtaining a better estimate for x^ and X 2 each time until the 
exact solution (or one within acceptable error limits) is obtained.* 
This concept can be extended to the general case by allowing AX 
and F to be n-dimensional vectors and J to be an n-by-n matrix. 

For the case at hand, the Newton-Raphson zero finding technique 
is applied to (3-8) . Except at the point of solution the results 
are not zero but 


P 


n 


= c 

n 


e 

n 


(3-17) 


* This assumes convergence will occur. In general, the proce- 
dure will converge to f(£) - 0 if f(x), f’(x), and f" (x) are 
continuous on some interval containing £, if f ' (&)£ 0, and if the 
starting point is "close enough" to the root £ . These conditions 
cannot be assured for a general nonlinear function or with randomly 
chosen starting values. See Conte and de Boor [17]. 
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The desired result, of course, is for p^ to he equal to zero. 

Note that the values are produced by the system of Figure 3.1 
operating on the- assumed input values of e. Thus, for two input 
variables x^ and x 2 > the desired form of (3-17) is 

P 1 = f^x^x-g) - x 1 = 0 (3-18) 

^ P 2 = f 2^ x l ,x 2^ " X 2 = ° 

where x^ and comprise the input signal e. 

k ^ 

1 * ♦ k ! 

The Newt on-Rap hson technique is applied to the variable p , 
resulting in 



P 


2 


+ [J] 



thus > 



The Jacobian matrix J is given by 


9p l 

9p l 

9x l 

9x 2 

9p 2 

9p 2 

8x i 

9x 2 


(3-19) 


(3-20) 


(3-21) 


The partial derivatives are computed by 



hk 


3P X 9 



and 


9p l 

9f l 

. ^ 

9x 2 " 

9x 2 

9x 2 





9x 2 


3P 2 

9X 1 

9 

9x^ 

C f 2 (x l 


3f 2 

9x 2 


9x l 

9x l 



9p 2 


^2 

^ X 2 

_ 

3 x n 



2 

2 




_ 

2 

1 


9 x 2 



The Jacobian matrix may thus be written 



(3-22a)' 


(3-22b) " 


(3r22c) 


,(3-22d) 


(3-23) 
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The system of Figure 3. 2 can be equivalently represented as a 

single block as shown in Figure 3»5. Let the input now be CN and the 

output YN. Let these correspond to the complex components of the 

Fourier series representation of the input and output signals. The 

component parts of CN are: a , the dc component; b, , the imagi- 

nary part of the fundamental component; a ^ and b^ , the real and 

imaginary parts of the second harmonic, etc. The a^ component is 

assumed to be zero so that the input signal corresponds to a sine 

* 

wave at the fundamental frequency. In addition, the fundamental 
frequency <a is considered an input. 

The system operates on these inputs to produce the output YN. 
The components of the output are the same as the input except that 
a nonzero a^ component is possible. The frequency m . is not 
changed by the system and thus is equal to the input io . 

The input CN component parts, through the second harmonic, are 

represented by a^ , b^ , , and b^ (remember that the real part 

of the fundamental component is zero and there is no imaginary part 
of the dc component) and the output YN components are , 

A^ , and . Equations 3-18 take the form 


*The actual reason fdr setting the real part of the fundamental 
component to zero is to establish a fixed time reference. 




4S 

On 
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P 1 A 0 a 0 (3-24) 

P 2 = A^ - 0 (since a^ = 0) 



The Jacobian matrix then becomes 


3A 

ar-1 

0 

3A 

0 

9A 
• 0 

3a 2 

3A 

0 

3b 2 

3 A 

. » 0 
3o> 

9A 1 

SA i 


3A 1 

^1 

3a 

0 

9b i 

3a^ 

3b 2 

3ta 

9B i 

8B i , 

3B 1 

9B 1 

9B 1 

3a 

0 

3b 

3a 2 

8b 2 

3(0 


3A 2 

3A, 

± 

9A 2 

^2 

3a 

0 

3bi 

3a 2 - 

9b 2 

3(0 

3B 2 

3B 2 

3B 2 

3B ? 

3B, 

3a 

0 

9b l 

8a 2 

3b 2 

3o) 


(3.25) 


Notice that because the component is always zero, there are no 


3a, 


terms. This eliminates one column of the Jacobian matrix. 


3 

3to 


since to is an 


There is, however, a column corresponding to 
input. The matrix thus remains square. 

The partial derivatives in the Jacobian matrix are computed by 


numerical means. This is accomplished by computing the ratio of 
the change in each output to a change in the input, e.g. 
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9A AA n 

o a 0 

9 a Aa n 

o 0 


(3^26) 


As long as the'ideltas in the fraction are. small, the approxima- 
tion is sufficiently accurate. These values are computed, one 
column at a time, by obtaining the change in the output values as 
one of the input components is altered slightly. That component is 
returned to its original value and another is altered to calculate 
another column of values. The procedure is repeated until all 
values are derived. 

Once all of the elements of the J matrix are obtained, the 
incremental changes in the input values are determined? from 


r ~ 

’ 

— — 

> 

o 


p i 

st i 


p 2 

da 2 

= [j r 1 

P 3 

Ab 2 


P 4 

Aw 


1 

in 

Qu 

1 


(3-27) 


where the p*s are as shown in Equation 3— 18. These delta values are 
then combined with the CN components to obtain the new values of CN 
to be used in the next iteration. 


3.6 Iteration Limits 


The numerical difference method of determining the partial 
derivative values of the Jacobian matrix is valid if the difference 



49 


values are sufficiently small. In order to insure its validity, 
the change in each output value is monitored to assure that the 
change is not too large. The criterion chosen is based on limiting 
the maximum change to about 2.5% of the value of the output variable. 
If any of the outputs change more than this limiting value, the 
incremental size of the input producing the excessive output change 
is reduced by a factor of 2 and the outputs are recomputed. This 
process is continued until all output changes fall within the 
established limits. 

Although the Newton-Raphson procedure rapidly Converges to a 
solution .in most cases, it is possible that the procedure will not 
converge at all. The automatic step size generation feature may 
result in step sizes which are too large to produce convergence. In 
order to circumvent this problem, the iteration step size is con- 
strained to be no more than 1.5 times the last successful iteration 
step size. Also, for each iteration a cost function is computed by 

I 

COST = l p. (3-28) 

i~l 1 

This cost is compared to the cost computed for the previous iteration; 
if the cost has not decreased, the iteration step size is reduced by 
a factor of five and the cost is recomputed. If four attempts are 
made without a cost reduction, the last reduced stepsizeis accepted. 


*It should be noted that by limiting the step size the Newton- 
Raphson procedure is really a gradient procedure. 
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Recognizing that numerical noise may become significant for 

—8 

computed values less than 10 , provision is made to terminate the 

process when iteration step sizes are consistently of this magnitude. 

_8 

Four consecutive iterations with step size less than 10 will 

terminate execution of the code. 

The procedure is also terminated when an output limiting func- 
—8 

tion is less than 10 . This function is found by 


YNOR 



(3-29) • 


“8 

When this limiting function is less than 10 , all output components 

are assumed to be zero. 

The limiting function is also used to identify another stop-- 
ping criterion. If each component of the output differs from the 
corresponding input component by an amount which is less than 0.1% 
of the limiting function of Equation 3-29 the procedure is terminated; 
a limit cycle has been successfully identified. 

Several example systems were considered using the procedure 
outlined above. The results which were obtained are presented in 


the next chapter. 



CHAPTER IV 


RESULTS OF APPLYING THE HBA 

The Harmonic Balance Algorithm (HBA) discussed in the previous 
chapter was applied to several example problems. These examples 
include single-path, multiple-nonlinearity systems and multi-path 
systems; they incorporate both memory and memoryless type non- 
linearities. The results obtained are summarized in Tables 4.1 
through 4.6. 

The system of Figure 4.1 has two identical nonlinearities 
but dissimilar linear elements. With point A as the input, ; 
randomly chosen initial values were selected; the initial value of 
id was set at 0.5 radians per second with a fundamental component 
amplitude of either 5 or 10 units. With an initial input of 5 units; 
the HBA terminated when the computed input increment consistently 
remained less than the perturbation for determining the J-matrix. 
However, where the input was initially set to 10 units, a limit 
cycle was identified. The frequency was determined as 1.015 radian/ 
second -with a fundamental amplitude of 18.16 units. As might be 
expected from the asymmetrical nature of the nonlinearities, a small 
dc component of- -0.867 units was also -determined. Second harmonic 
cosine and sine components were determined to be -0.124 and -0.110 
respectively, while the third harmonic components were -0.819 and 
-0.374. (Only components through the 3rd harmonic were calculated.) 

With the order of the linear elements reversed by choosing 
point B as the input, identification of a limit cycle proved to be 



INITIAL VALUES 

Number of Harmonics = 
e(t) = 10 sin 0.5t 

Number of Harmonics = 
e(t) = 5 sin 0.5t 

Number of Harmonics = 
e(t) = 5 sin 0.5t 


Table 4. 


3 


3 


COMPUTER TIME 
(MIN: SEC) 

I 


HBA RESULTS 


1:29.9 e(t) = -0.867 + 18.16 sin.l-._Q15 t 

-0.124 cos 2.029t - 0.110 sin 2.029t 
-0.819 cos 3. 045t - 0.374 sin 3.045t 

2:04.4 Terminated On Step Size 


1 


:14. 7 


Terminated on Step Size 


HBA Results For System of Figure 4.1, Input At A . 


Ur 

ho 




Number of Harmonics =1 :22.4 

e(t) = 0.1 sin 0.5t 

Number of Harmonics = 1 :24.7 

e(t) =-0-. 1+ 2 sin t 


Table 4.2: HBA Results For System 


HBA RESULTS • 
Terminated on Step Size 


Terminated on Step Size 


Terminated on Step Size 


e(t) = -0.064 + 
1.811 sin 1. 054t 


Figure 4.1, Input At B . 


Ut 

w 



INITIAL 


VALUES 


COMPUTER TIME 
(MIN: SEC) ' 


HBA RESULTS 


Number of Harmonics = 1 !|04.6 

e(t) = 0.1 sin 2t 

Number of Harmonics =1 :18.9 

e(t) “ 4.0 sin 2t 

Number of Harmonics =1 :13.1 

e(t) - 0.5 sin 3t 

Number of Harmonics =1 :15.0 

e(t) = 10.0 sin 8t 


Amplitude = 0 " 


e(t) = 5.43 sin 4.00 t 


e(t) = 0.410 sin 4.00 t 


Amplitude « 0 


Number of Harmonics = 3 

:53. 0 

e(t) - 

5.33 sin L 

u09t 



- 

0.026 

COS 

8.19 t 

e(t) = 4.0 sin 5t 



0.024 

sin 

8. 19 t 


+ 

0.488 

cos 

12.28 t 



+ 

0.239 

sin 

12.28 t 


Table 4.3: HBA Results For System of ' Figure 4.2 With NL2a . 



INITIAL 


VALUES 


COMPUTER TIME 
(MIN: SEC) 


HBA RESULTS 


Number of Harmonics = 3 
e(t) - 1.0 sin 2. 5 t 


:52. 1 


e(t) = 0.502 sin 3.58 t 

+ 0.12xl0” 4 cos 7.15 t 

- 0,98x10“5 sin 7.15 t 
+ 0.0336 cos 10.72 t 

- 0.0218 sin 10.72t 


Table 4.4: HBA Results for System of Figure 4.2 with NL2b 



INITIAL VALUES 

Number of Harmonics = 1 
e(t) = 0.4 sin 2.5 t 

Number of Harmonics - 3 1:15.7 

e(t) = 0.4 sin 2.5 t 



Table 4.5: HBA Results For System 


HBA RESULTS 

e(t) => 0.467 sin 3.91 1 


e(t) ® 0.454 sin 4.04 t 

- 0.5xl0 -5 cos 8.07 t 
+ 0.1x10”^ sin .8.07 t 
+ 0.0228 cos 12.11 t 

- 0.0230 sin 12.11 t 


Figure 4.3. 


Ln 

ON 



INITIAL VALUES 


COMPUTER TIME 
' (MIN: SEC) 


Number of Harmonics =1 :28.0 

e(t) = 0.25 sin 0.11 t 

Number of Harmonics =3 2:08.1 

e(t) = 0.25 sin 0.1 t 


Table 4.6: HBA Results For System 


HBA RESULTS 


e(t) = 254.7 sin 0.100 t 


e(t) = 265.1 sin 0.0932 t 
+ 0.00166 cos 0.186 t 

- 0.00286 sin 0.186 1- 
+ 22.73 cos 0.280 t 

- 13.84 sin 0.280 t .. 


Figure 4.4. 


Ln 
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more difficult with randomly selected values. A number of frequency 
and fundamental amplitude combinations resulted in termination 
based on the step size criterion. When a small dc component was 
added, however, a limit cycle was identified. Selecting the input 
as to = 1.0 , DC = -0.1 , and fundamental = 2.0 resulted in conver- 
gence to to = 1.054 , DC = -0.064, and fundamental = 1;811. (Wo har- 
monics were calculated in this example. ) Since these two closed 
loop systems are identical except for the point chosen as the input, 
they should have identical limit cycle frequencies. Notice that the 
frequencies determined by the HBA agree Within 5%. 

Davison and Constantinescu (14] considered several examples in 
illustrating their method of limit cycle identification. These 
examples were also examined using the HBA. Limit cycles were iden- 
tified in each case. 

The systems of Figure 4.2 involves three nonlinearities in a 
single path; with NL2a in the system, several combinations of ini- 
tial to and fundamental amplitude were selected, ranging from to = 
2.0 to 8.0 and amplitude from 0.1 to 10.0 . When considering only 
the fundamental component (i.e. , -no harmonic components) limit 
cycles were identified at to = 4.00 rad/sec. Two limit cycle ampli- 
tudes were possible, 0.410 and 5.43. In general, for to t 6.0 , or 
with amplitude less than or equal to 0.5. and to = 2.0, the HBA 
terminated with the amplitude equal to zero. This is, of course, a 
possible solution although it yields no useful information. 

Limit cycle identification was also achieved when three 
harmonics were considered. An initial input of 5 rad/sec with a 
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fundamental amplitude of 4.0 produced a limit cycle frequency of 
4.01 rad/sec. with a fundamental amplitude of 5.33 units. Second 
harmonic components were less than 1% of the fundamental component, 
but the. third harmonic components were somewhat . larger . A third 
harmonic cosine component of 0.488 and a sine component of 0.239 
were identified. 

The Davison and Constantinescu.: method identified a limit cycle 
at 4.00 rad/sec. with a fundamental amplitude of 0.58 or 5.32. 

While the frequency agrees with that determined by the HBA, the 
amplitudes are somewhat different. This may be due to inherent 
inaccuracies of graphical techniques as employed by Davison and 
Constantinescu. 

The second example considered by Davison and Constantinescu 
replaces the second nonlinearity by a relay with hysteresis-. (The 
hysteresis gap is + 0.2 and the amplitude is + 1.5.) A limit cycle 
with an amplitude of 0.50 units at 3.61 rad/sec. was identified. 

The HBA, considering three harmonic components, identified a limit 
cycle of 3.58 rad/sec. and a fundamental amplitude of 0.502 units. 
The second harmonic components were less than 0.01% of the funda- 
mental, but the third harmonic components were larger; the third 
harmonic cosine amplitude was 0.0336 and the sine amplitude was 
-0.0218. 

A modification of this example was also considered using the 
HBA. The first nonlinearity of Figure 4.2 was altered slightly in 
addition to the replacement of the second nonlinearity by the relay 
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characteristics. The system of Figure 4.3 resulted. This pro- 
duced a slightly different limit cycle, one with a frequency of 
3.91 rad/sec and a fundamental amplitude of 0.467 units. When 
three harmonic components were considered, the limit cycle fre- 
quency shifted to 5.04 rad/sec. The resulting fundamental 
amplitude was 0.454 units. The second harmonic amplitudes were 
again negligible, but the third harmonic components were not, 
having amplitudes of 0.0228 and -0.0230 for the cosine and sine 
values respectively. 

The last example considered was a parallel path system as 
shown in Figure 4.4. The HBA was applied for three harmonics, an 
initial frequency of 0.1 rad/sec., and a fundamental input compo- 
nent of 0.25 units. A limit cycle was identified at 0.0932 
rad/sec. with a fundamental amplitude of 265.1 units. Second 
harmonic amplitudes were less than one unit, and third harmonic 
cosine and sine values were 22.73 and -13.84 respectively. Davison 
and Constantinescu identified the limit cycle at 0.1 rad/sec. with 
an amplitude of 254.8 units. When the HBA was limited to a single 
harmonic component, the limit cycle was determined at &> = 0.100 
with an amplitude of 254.7 units. 

The agreement of the HBA results and the Davison-Constantinescu 
results support the validity of , the HBA metho.d, but the examples 
also show that harmonics can alter the limit cycle. In each case 
where multiple harmonics were considered, both the frequency and 
the fundamental amplitude were somewhat different from the single 


harmonic case. 





















CHAPTER V 


CONCLUSIONS AND RECOMMENDATIONS 

Techniques for determining limit cycles in nonlinear systems 
rely heavily upon describing fuction analyses which considers only 
the fundamental component of the Fourier series- representation of 
the nonlinear element output. Thus any contribution that may be 
made by other harmonic components is. ignored. No practical method^ 
of analysis considered multiple harmonic components in systems 
containing several nonlinearities or forward paths. Since the 
effect of these components is not always negligible, a method 
of analysis was needed which included them. In an effort to meet 
this need, the Harmonic Balance Algorithm was developed. 

The HBA can be used to identify limit cycles in multiple non- 
linearity, multiple path systems. It cannot, however, be generally 
applied in a random manner. Random selection of initial input v i 
values may result in the identification of a zero amplitude limit 
cycle (a trivial case) or termination of the iterative process due 
to insufficient step sizes. Thus it appears that the initial input 
values should be intelligently selected, perhaps based upon the 
results of some other method (the Davison-Constantinescu techniques 
being the most general) . In the event that a limit cycle is not 
identified by the describing function or other techniques, but 
indications are that one almost exists (i.e., the curves almost 
intersect) , the HBA would be useful in determining if a limit cycle 
exists when harmonic components are present. 
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At any rate, the HBA cannot be profitably used by an individual 
unfamiliar with nonlinear system analysis. Its success depends in 
large measure upon the user T s insight and ability to select the 
initial input values; it should not be considered a general purpose 
tool to be used to circumvent analysis using other available methods. 

While the HBA sucessfully identifies limit cycles, its success 
is highly dependent upon the selection of the initial input values. 
While methods exist for determining possible limit cycles when only 
one or two harmonic components are considered, they are generally 
very tedious and somewhat error prone. Techniques for simplifying 
these methods, as well as new techniques, would be appropriate. 

One specific area which. could be examined is the Davison- 
Constant inescu method. A computer algorithm for generating the 
required curves or the data for plotting them would substantially 
enhance the attractiveness of the method. The limit cycle or near 
limit cycle values suggested by this method could then be used as 
starting values for the HBA in examining the effect of harmonics. 

Further work may also prove fruitful in the area of automatic 
step size generation. In its present form the HBA may converge to a 
nonzero minimum of the input-output difference function. When 
this occurs, the algorithm terminates because further refinement 
of the input values will not produce the desired results. Some 
method of automatic input selection may be possible for these cases 
which would extricate the algorithm from these local minimum traps, 
allowing the zero finding process to continue. 
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Using the HBA as a base, additional coding may allow 
consideration of a more general system. In its present form the 
HBA does not allow cross-coupling between forward paths, i.e. 
feeding the output of an element in one path to the input of an 
element in another path. Also, no provision is made for internal 
feedback loops within a path. Both of these variations appear 
possible using the general concept upon which the HBA is based. 
Further efforts in this area should prove fruitful. 



.APPENDICES 
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APPENDIX A 

HAEMONIC BALANCE ALGORITHM USER'S GUIDE 

The Harmonic Balance Algorithm (HBA) is intended to be a 
nonlinear system analysis tool. Its success depends upon the 
intelligent selection of input data; an individual untrained in 
nonlinear system theory is unlikely to be able to use it beneficially. 

With one exception, data is entered in free format form; i.e. , 
data entries are not sensitive to columnar location on the input 
data card. Multiple data entries on a single line or card are 
separated by a comma, denoting the end of one data field and the 
beginning of another. 

The single exception is the first data card, which designates 
a number of iteration variables and limits. Blanks or zeroes in 
the fields defining these inputs will cause predetermined values to 
be selected. These values are defined in Table A.l. 

The significance of input variable NFLAG requires some further 
explanation. Under normal use only a limited amount of output data 

I * 

is required. As shown in Appendix C, this includes the iteration 
number, the size of the largest pertubation used’ in calculating the 
partial derivatives in the J-matrix, the cost figures, the maximum 
input increment, and the new values of input to be used on the next 
iteration. The user may choose to output additional data by select- 
ing larger values for NFLAG. Each larger number results in more 
intermediate data being printed, such as the J-matrix values, the 
outputs of the nonlinearities, the linear system gains, etc. When 

precomns page blank-not & 
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VARIABLE FORMAT COLUMN DEFINITION 


DELTA 

E 

1-10 

Perturbation used in calculating the 
J-Matrix (Default = 10“®) 

NPMX 

I 

11-20 

Number of points used to describe the 
inputs function (Default = 301) 

ITEKNO 

I 

21-30 

Maximum number of iterations in calcu- 
lating the output of nonlinear element 
(Default = 100) 

EPS 

E 

31-40 

Secant method stopping limit (Default 

= 10“5) 

EE 

E 

41-50 

Minimum step size limit for termination 
(Default - 10-8) 

SIZE 

E 

51-60 

Initial Maximum input increment 
(Default « 0.1) 

NFLAG 

I 

61-65 

Output -Data selection (Default « 0) 


Table A.l. Input Variables on First Data Card. 
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NFLAG is set to 3, all WRITE statements are used; intermediate 
calculations and subroutine entries and exits are indicated. The 
latter case may prove useful in helping one to understand the. 
intricacies of the program. 

The remaining input data can best be explained by use of an 
example. Each input variable is defined in Table A. 2. and will be 
identified as the sample input data in Table A. 3 is examined. Note 
that each input data card or line is numbered. These line numbers 
are not actually present on the input data card , but are included 
here for reference purposes. A flow chart depicting the data input 
process is shown in Figure A.-l. An -explanation of each line follows. 

Card Is DELTA, NPMX, ITERNO,. EPS,. EE, SIZE, NFLAG. This is the 

only input data not using free format. A blank or a zero 
for any value results in default selection. 

Card 2: Number of .harmonics (NOHAR) and number of parallel paths 

(HP). 

Card 3: Initial guess at limit cycle frequency (WO) . 

Card 4: Initial guess at limit cycle amplitudes (HAR, entered in 
order of DC, fundamental, second harmonic cosine, second 
harmonic sine, third, harmonic cosine, etc.),. 

Card 5: Number of nonlinearities (NNL) in each of the paths (path 

one, path two,, etc.). 

Card 6: Number of breakpoints (NBKPTS) for the nonlinearity and 

the nature of its symmetry (NLSYM) . The number of break- 
points is one more than' the actual number since an 
artificial point is inserted at the far left of the 
nonlinearity characteristic. 

Card 7: JD, the nonlinearity curve number. (See Figure A. 2.) 

This is always 1 for memoryless nonlinearities. If the 
nonlinearity is a memory type JD is 1 or 2 depending upon 
which curve is being examined. 

Card 8: Abscissa coordinates for the nonlinearity breakpoints 

(EBKPT) beginning at the artificial left coordinate and 
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VARIABLE 


DEFINITION 


NOHAR 

1 NP 
. WO 
HAR 
NNL 
NBKPTS 

NLSYM 

JD 

EBKPT 

DBKPT 

PWGAIN 

NH 

AB 

MORDER 
NORDER 
AA . 

BB 


*of the form 


Number of harmonic terms' to Be considered in addition 
to DC term 

Number of parallel paths 

Initial guess at limit cycle radian frequency 

Initial guess at input fourier series components 

Number of .nonlinearities in each path 

The number of breakpoints in the piecewise linear 
representation of the nonlinearity 

Nonlinearity symmetry characteristic (0 = asymmetri- 
cal, 1 = symmetrical) 

■v , \ 

Nonlinearity curve designation (used to identify 
curves on hysteresis type nonlinearities) 

Nonlinearity breakpoint abscissa 

Nonlinearity breakpoint ordinate 

Nonlinearity piecewise-linear slope 

Number of common breakpoints on hysteresis type 
nonlinearity 

Hysteresis nonlinearity common breakpoint abscissa 

Order of numerator of linear transfer function 

Order of denominator of linear transfer function 

Coefficients of linear transfer function numerator 
terms* 

Coefficients of linear transfer function denominator 
terms* 

V” 4- + A 3 S n - 2 + ■ . ■ + A n S + A^ ■ . 

+ BjS 1 "- 1 + B/* 2 +....+ B m S + B^ 


Table A. 2: Input Variables Using Free Format 
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Card or 
Line Number 


Data 


1 

(Card 1 

was left blank in 

this exai 

2 

3, 

0 



3 

2.5 




4 

0.0, 

1.0, 

0.0, 

0.0, 

5 

3 




6 

3, 

1 



7 

1 




8 

-l.E+05, 

-2.0, 

2.0 


9 

-3.0, 

-3.0, 

3.0 


10 

0.0, 

0.0, 

1.5, 

0.0 

11 

0 




12 

0 




13 

1, 

2 

, 


14 

5.0, 

20.0 



15 

1.0, 

3.0, 

5.0 


16 

4, 

1 ' 



17 

1 




18 

-l.E+05, 

-0.2, 

0.2, 

0.2 

19 

-1.5, 5 

-1.5, 

-1.5, 

1.5 

20 

0.0, 

0.0, 

0.0, 

l.E+20 

21 

2 




22 

-l.E+05, 

-0.2, 

-0.2, 

0.2 

23 

-1.5, 

-1.5, 

1.5, 

1.5 

24 

0.0, 

0.0, 

l.E+20, 

0.0, 

25 

0 




26 

2 




27 

-0.2, 

0.2 



28 

2, 

3 



29 

1.0, 

3.0, 

5.0 


30 

1.0, 

7.0, 

16.0, 

12.0 

31 

3, 

1 



32 

1 




33 

-l.E+05, 

-3.0, 

3.0 


34 

-2.5, 

-2.5, 

2.5 


35 

0.0, 

0.0, 

0.83333, 

0.0 

36 

0 




37 

0 




38 

o. 

1 



39 

10.0 




40 

1.0, 

4.0 




0 . 0 , 0.0 


0.0 


Table A. 3: Sample Input Data for Single-Path, Three 

Nonlinearity System With Hysteresis 
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NREAD. DELTA, NPMX, ITERNO, EPS, EE, .SIZE, NFLAG 

~n ' ' 

READ NOHAR, NP 


DELTA 

NPMX, ITERNO, EPS 
„ S>> ^EE, SIZE 


Set Zero 

Value to Default 
Value 
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progressing to the. right. The artificial breakpoint is 
set far enough to the left so that it is never crossed. 

Card 9: Ordinate coordinates for the nonlinearity (UBKPT), begin- 

ning at the artificial left breakpoint and proceeding to 
the right. 

Card 10: Slope of the nonlinearity (PWGAIN) to the left of the 
associated breakpoint. - There is always one more slope 
than breakpoint, the last being the slope to the right of 
the rightmost breakpoint. The slope to the left of the 
artificial breakpoint is set to zero. Vertical slopes 
are indicated by a very large number (i.e., 10 ^ 0 )^ 

Card 11: JD. If JD = 2, cards 8-10 are repeated for curve 2 of 

the hysteresis characteristic. JD = 0 indicates the non- 
linearity has been fully described. 

NH, the number of hysteresis curve common point. NH = 0 
for a memoryless type nonlinearity. 

i 

Order of numerator (MORDER) and denominator (NORDER) of 
the linear transfer function. 

Card 14: Coefficients of the numerator polynomial (AA) of the linear 
transfer function, with the coefficient of the highest order 
term first. 

Card 15: Coefficients of the denominator polynomial (BB) of the lin- 

ear transfer function, highest order term first. 

This completes the entry of data for the first nonlinear-linear 

pair for path 1. If other pairs remain in path 1, the process 

(cards 6 through 15) is repeated until all path 1 data is entered. 

Data is then entered for each successive path until all data is 

entered. 

Cards 16 - 20: Same as cards 6 - 10, but for nonlinearity 2, 

curve 1. 

Card 21: JD = 2, thus breakpoint information for curve 2 

of the hysteresis nonlinearity follows. (See Fig- 
ure A. 2) 

Cards 22 - 24: Same as cards 8-10, but for nonlinearity 2, 

curve 2. 


Card 12: 
Card 13 : 
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Card 25: 

Card 26: 
Card 27: 

Cards 28 
Cards 31 
Card 36: 

Card 37: 
Cards 38 


JD ?= 0, so all breakpoint information for this non- 
linearity has been entered. 

Kumber of hysteresis curve common points (2) . 

Abscissa common breakpoints for the hysteresis 
curves (-0.2, 0.2). 

30: Same as cards 13 - 15, but .for second, linear element, 

35: Same as cards 6 - 10, but for third nonlinearity. 

JD *f 0, so all breakpoint information for third 
nonlinearity has been entered. 

NH = 0, so no hysteresis common points exist. 

40: Same as cards 13 - 15, but for third linear element. 


This completes the data entry for path 1, and: since only oneo 
path is used in this example, all data has been entered. For multi- 
path examples, data entry for path 2 would follow card 40 and would 
follow the pattern begun at card 6. 

Once the data entry phase is completed,; no further user input 
is required. The HBA outputs identify possible limit cycles as 
determined from the initial inputs. Note that failure to identify 
a nontrivial limit cycle does not mean that no limit cycle exists; 
it simply means that, if one does exist, it has not yet been 


identified 
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APPENDIX B 


FORTRAN CODED LISTING OF HBA 
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D£LT<ILIM) 

?70 

CONTINUE 



DO 2E0 1=1. I L IM 


:-:£u 

DELTCL ’=1.51 BEL 

fiO 


crfcf j mhipi:-: 

!“0P £EPQ POMS 


C IlfFHTIi-iE? WHICH ROMs mRE ?£R0 BY ii 0 

H.M. l.’niKl. l)-l. 

00 P'UT k=ILIM. 3 -t 

k- 1 J--i r J'.K k . 

TOO C UN I HU IF. 

BO 3tO 1H=1> ILIM 
SVi iL=C. 0 

nn 300 JM“Ij ILIM 

UVhL-BVi A rtJi IM. JH >v 

300 com IhUE 
•JFLAG 1 ' IM)“1 

IF <'I!'1.£0.2> OJ TO 310 
, IF * BvY)L« ! " 1 . 1 . fl -1 0 ' 

310 CONTI HUE 

IF • NFLrtG, GE.'t J WRITE <6,770' 

BO £30 K=l- TI.IM 

IF \NFLhG. GE. J i UR I IE <.6. 66 U j h. * 

330 IF <NFl rtG.GE. J i WRITE <6. 650) <«JtK. JV’>, JV=1. ILIM) 

IF CHFL.uG.uF. 1’’ WRITE <6.320’ <BF.lT<I>. I~i. Il.IMv 
C ELTMINhTF ROWS h!1D COLUMNS CORRESPONDING to HERO f.Mji.t SUOSCRI 

110 350 n-ILIM. 1.-1 

IF ■•jrLHC.'kl U-.C.I.' uo.lu 35 U 
BO 330 I 2«1 .11 li'i 

BO 330 r.3;jFl..> ILIM ' 

hJ* 1 2- 1 'J'T'.-B F3+1 1 , 

3^0 CONTINUE 

BO 340 K3-1. IUM 
I'U 340 K?"Fi ■ ILIH 

nB.f'2.k3T. TiJ-'l 2i-l. H; 

3.5 0 CONTI 1 " IF ‘ * 

■'*>0 Cj w 1 i NuF 
\ BIG=0 

DO 360 lBUMP'- 1 . ii.rri 

C LtJiU IS THE hUnCtR OF NON-FkRu RO'-'S 

F£-:IG=k>-: IG-i Jf Li iC. IFUnP.' 

360 CONTI NOR 

1 . WRITE hujUSTB.' J MhTkJ.4 

IF (.HFLisG.GE. 1 1 WRITE <6.600’' 

DO 3 c Q IS- 3.. Kb Ik 

IF • nFI.hG.GR.i ■ MR I '6 '‘6.66'.;' Iu 

IF ' NrLOC.Q-. i WRITE 65 O'. m«J< 111. JVI .. .VI* j • KBI'.F* 

370 CONTINUE 

r n:.u inotrsc u* j MuTPT : 

ChLL IthT/NY 'iTJ.J-FiTG hJINV. IFR. 25 ■ 

IF ''NFLtU,.GE.3> WRITE <.6.610> 

IF (.IER.EO.,? ' GO fU_ 460 
I \- • Nt-LnCk GE. 1 ’• WRIT E ■' 6 ■ ?6 0 ’• 

DO 380 I=i.»I r . 

IF (-NT-1 hG.GE.O WflfE ih. 67Cv 1 

IT ‘Tl' LnG. GE. 1 ’> WRITE ''S- cOC < i/iJiHVt L- JL/. JU-3 . Kbit ■ 

380 CONI I NOR 

r Cm! CUuT 1 1' U'UC I I Oil V FOR US! IN UNOtHG li 'l'll! BE I inS 

Cpu. Fi'iii • -fOr : 1 v. w.Ir-i. tnR. Ji cmG. Chi . jcbi-.lI'* ''2 . 1 ■ 

IF CNFL: 3’* WRITE f '6. S2Q ■ 

USING ChLCl.I fiTEn RUNlTiOM rtMD J irlv'EuSi'. BE 1 1 : Ml N:-. BE- in 
1 .. vulUEt i-'OR Till-. INPUT 

ChLL DEI .* ('ii-LCU, OJTNV. ILIM. FV- RN'G. Jl-1 h£- D. SI 
IF ■ NT* Li :G.Gi-.3 ■ <.6-6300 
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DEI ERMINE Mi. *'l VALUES UF Crl FUK USE TH ITFRhTIOM 
ICGST=0 
C03T1=CGST2 
Ml *1.10 

DO 390 L-l,JrtAX 
390 • GNULXNJ'.L / 

GO TO 420 
400 WRITE ^?80.« 

IF acosr.GE.-P GO TO 40 
DO 4l 0 K“ij ILtM 
• s i 0 BFLCtKK'-BELEN'L) '5. 

D"D^5. 

KChT=KCHi-i 
WRITE <6* ?90> KENT 
IC08T=ICGST+1 
420 MO=Mi-DELCN' ILTm 
WRITE MO 

CHKn=CNia:>-GMFLX^OFLCNa>, 0. 0) 

DC-REAL* Crll e J )' 

WRITS: <t~* 73 U’ 1 DC 

CMK2?«cm <2 t-CMPLX^ 0.0. BELCH f 2 * 

' hUND“-2 . *hIHAG<CHK2) ) 

WRITE ce. 740; FUND 
JF UMhX-LE- 2.) GO TO 440 
DO 430 KJ-3- JrlwX 
KR=<2 *KJj-3 
KI-<2+FJ>~2 

CrJKKJj-CHl'.P J '-CItPlX^DELCN' KPi. BELCH <KI>> 
HARC=2. +REOiy CHI < K J> T 
H ARS=~2 .Uni MAG c C f II C K J ^ 

JK-kJ-I 

WRITE u5,750> JK,HAPC,HHRS 
430 CONTINUE 
440 kCHT-KCm f+i 

IF c KCN T . GT. 1 00') STup 
YMAN- ABo'DELT • 1 * / 
nu 450 K=i* II IM 

450 YMi iX* AMAX t t YMhX . A&S ■* DE L T < J > A 

IF Ol.LT.YMfO NB-ML+1 
Ih UJ-Gt -YMA>.> rtfNQ 
IF 'ND.GT.2' GO TO 480 
3 I 3F- 2,^0 

I h \Sl2£»L I . ££j nSD»HSC+ 1 
IF » 0I2E,GT.EE • HSOO 
IF * IISU. GT. 4>‘ GO fO 470 
if'FITF 1 r.wlO* KCHT 
WRITE *'^8lsj* 

GO TO 30 

4 *S0 WHITE '.C> S0 T V' 

STOP 

470 URITh <F.S30« 

STOP 

480 WPITF: <6.£4 0> 

STOP 



490 n iR ITE *' 6 . *7 m « Ti«“ 

WRITE <6.880' FUND 

IF « JMhH.LE.2^ G'j TO 510 

DO 500 J.K3- JMAX 

Hj iR02 . * PEAL <CNKJJ> > 

hmrs=«2. *himagu:iu <jj/> 

son WRITE <6. 9nn> JH, HARE* HuPS 
510 MPIiE *6. 89 1U Mu 
STOP 

52 U WRITE ?I0^ 


530 KORilw T >’2 0>i, 'PtTniER MAIN FROM RMAIT> 
540* FORMAT <2 OX- ‘REG ITER MAIN FROM OUTPUT R‘ 
550 FORMAT i2u::. 'REENTER MmIN FROM FUHi"; 
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M 0 FORMAT '2X, 'CNI- M3>0= J ,2<Gli.5,5X>/'2X. 'CMD'M3» '-2<G11.5/ 
'DELTiMA ' )= ',G11.5,5X, 'WO= ',G11.5> 

570 FOPMAT <20X. ' REENTER MAIN FROM OUTPUT V, REAL VAL I NCR " 1 
580 FORMhT -.5X- 8>‘Gi3 .5, 3X>) 

590 FOPMhT 120X- ' REENTER MAIN FROM OUTPUT V, IMAG VAL INCRO 
600 FORMAT C20N. ’.RECHTEP MAIN FROM OUTPUT V, W INCREMENTED-’ ') 

610 FORMhT C20X, ’ REENTER MAIN FROM HaTINV' 

620 FORMAT - SOM, * REENTER MAIN FROM FUN’i” > - 
630 FORMhT -'SON. “ REENTEP MAIN FROM BELC’ > 

640 FORMhT <-3 CM, '■**■* STEP SIZE SMhLlER THhN PARTIAL INCREMENT ■»**') 

650 FORMAT < ’+' . 10M. 8-.E1 1 . 5, 2X/> 

660 FORMhT *'1M- * J< M2- ' . *> 

670 FORMhT <2X. * JINVC J » 12- ■’->= •’> 

680 FOPMhT •■ION. ' JINV DOES NOT EXIST.' } 

690 FORMhT CviON. 'hD JUSTED J MATRIX’, ✓> 

700 FORMAT ■:■+’. J5M. 87051.5. 2X>> 

710 FOPMhT 40X. 'tt* F'UH NUMBER', 13, * FOLLOWS UP) 

"720 FOPMhT ' 'V40X. -'NLM MO = ', &1S.7J 

7o0 hOPMui •:^40X, • NEW LC COnPONEHl = ’ .G12.7> 

740 FORMAT 4 OX, 'NEW FUNBAMENThL SINE COMPONENT = ' <0,12.7 > 

750 FORMAT ■•V4UX. 'HEM COMPONENTS FOP HARMONIC NO. M2. COSINE - MI, 
li2. r - 5X- SI At - •' . Gl2. 

760 FOPMhT 

7 70 FORMAT <30X. 1 J MATRIX’’ > 

“t’O FOPriiT 1 OX- 1 1 1 r * Fnli Eli COS! TtST *n u7j/3 

790 FGFfviT kUX’.- ‘ KFri hT ChLCULhT 10! I FROM PUN IS'- 

800 FORMAT -.ROM. '4*4 COST 1 - '.GJ 2.7. ION. ’ COST 2 - ’,G12.7. ■{.»’) 

810 FORMAT -40X. ’ Nt'il MhX. PaRIIhL FERTURBhTXuN - J • Gl2. 7.- 
S20 FOPMA’I UN. 1 NEXT Pi.RTIhL INCREMENT =" 8-.E1 1.5. 2\0> 


20X- ■ H DEN. Ta STEP SIZE CONSISTENTLY LESS TNnN iO-8 

jo.-.. coh ,s Lt.. cm*: ■ 12 , r ■ - * .5 <gi?, fjtvr<> 

4 OX. ’PI -N NUMBER 1 FOL LOMS 1 ) 

4 OX. ’MAN. PiiPHAL PERTUPBhI ION = * ■ G12.7 ■ 
no - 1 1 Gi2. 

• - '40.- •’ FUNDhMFHI ill. SINE COMPONIIN T 
890 FORMAT >! - 40X. ’LIMIT C i'CLF FPL OUENC.Y - ‘Gl?.. 
fn‘n'1 E»" iKf]' :T \ ■ i'iTP|pn-'iEM ! c FOP a -•Mm-Tm T C ijil, ■ 

?;§X. STNE = '’. Gjl.7/' 

910 FOPMhT <>, 30X- ' PcTUPN EP.%OR FROM FUMY’ / 


830 FOPfinf 
84 0 FORrlH ; 

850 FORMAT 
060 FOPi’ih'i 
8 ? 0 t ORHA ! '■ rr 
880 FOPMhT 


’ < G 1 . 7 > 

FhDIANS SEC 
COSINE “ 


Gi2. 


KHO 

rUNCTTON E • r- 

complex ch-:u.-, pcii 

COMMON mIOPT' ON -OMN' 1.10 ✓PQMH • UMAX 
E-REhL (.CfK 1 J 
DO 10 I --2, JMhX 

PCM-CN- I,-tCEXP<CMPLX<U,.FL0Ar-'ji-ljU.)0'»r'7 
1 0 E-E+RFAL-'Pi.N KONJi'i- PCN ' I 
• : E 1 URN 


END 


subkoi M inr oc- ! i ‘U t t-ei i H t s i he - iu rpu ! or a php all c l 

SYSTEM ’.SHU N PhPm! LF’I. Ph 010 AND M NONLINEAR 1 1 J f -S. 

subpo 1 - r • ,j t i.i-j r ru r < yo-jt ■ rut . 

COMPLEX CM- n r ,C. 1 x ") Ou !\ U 1 . TN 13 1 
COMMON 'OPT'- NP, NNL-5.- .'MOPT,' CM .'POMN' JMkX 
COMMON NFLhG 

IF -nFlhG.GChS - v’ilt i00> 

DO 10 J=1.JMhN 

IF -NF! hG.GF. 3J UPITE u5, lJC) i.NK<JJ 

i o S’o-j 1 1 . i.-=c Mr: n. o, o. o,' 

i ;pp- t 

20 DO ot 1 J 1. Ji-AX 
30 CN- J> !S CNt >’ A 
NIT-1 

40 CONTINUE 

ChlL OUT *.iH, ni 1 ■ NPF‘ - 
IF -.NFLhG Ch.3- I'^IU- 100/ 

NIT= NIT i-J 
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If- C N IT . 6T . NHL. < NPP ) ’> GO TU 60 , 

IF tNFLwG.GE.2) WRITE <6> 140 j NIT> NPP 
DO 50 J=l. JM hX 

IF • NFl ftG.GE.S) I.IPI1E <6, 130'- YN<J> 

50 CIK-O-YN 1 J> 

GO TO 40 
60 NPF-NPP+i 

do ro j-ljmh:': 

?0 YOUT <J''=YQUT < JVi-YHfJ' 1 

IF 'MPP.LE, HP) GO TO 40 
DO yO J=J . JMPX 
80 • YOUTCJ.^t.-n+YOUn Jj 
DO 90 J=i, JriOw 

90 IF <HPl hG.GE.1) I-IP1TE <$. 1500 J. YOUT <J> 

IF ■ HFLPG.GE.3; l.'PITF <6. UO^ 

PETUPN 

100 I-ORMhT -RON. •' ENTER OUTPUT < 
iiu FC'Pfif'. ! >.2 OX.. ‘ l.EmVE uinPOT' 1 
120 FQPhAT ''SON'. 'PEEN1EP OUTPUT) 

130 FOPilnT C' 

< 4 fi FOPf>k!i 1 2 OX/ ’ COnPI-EX Iti p UT FOR NONLIKEhRITY^ 13* 
150 FQPf'iOT ' 1 0>:.. 'COMPLEX YOUIYMS, J >“'/ 2GJ5.7) 

End 


iO 

20 

30 

40 

50 

60 

70 

SO 


90 

100 


SUBROUTINE ORBFP • TO. JJHBNB. PIMDMT '• I I ; 

0 INFUSION KCnT 0500>. KDOGcbOu,'. IfcivbUU).. J 1 ND..B^. 

in mansion tsn^.soo:* 

COMMON NFLhG . , n . 

IF UJFLhG.&E.3> WRITE '.6.90; 

DO £0 J=i • II 
DUMB- 0 . 0 
DO 10 1=1,11 

IF '.DMflB.Gl .TB<D) GO lO XU 
BUMB-TBU.' 

f=i_ 

CON I /NUt 
K DOC f J --K 
TCl'Ts < ---TB«.L'j 
DO 30 1=1.. II 
F=T DDG< T ' 

TB:':<n+i-i''=-TBCio 
DO 40 I4JI 
TB''I.'=TBXa) 

no so i--i, ii 

\ - F TinG ■ i :■ 

KCiTT I I + 1- I ) =J I NDXB ( K > 

DO 60 Tl.- II 

JIMDS’E? 1 I ‘-KChI I 
DO 70 1=1, II 
I- =7 DOC ( I> 

l.cmT c i r n -i '=fihd>:t<ki 

DO 80 1=3. II 

KIMBXT > T=KNT*I> 

IF ‘HELhf..GF..3> I.PITE <6.100- 

PETUPN 


fop.tht bso::. 'enter order > 

FOPF'hT '-20'.. 'LEuv't ORDER ’ j 


END 


Pl-tlH* j 13) 


i'Cu.n K'INDXT t500> 


C 
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SUBROUTINE OUT <.YN, N. NPP) 

SUBROUTINE OUT COMPUTES THE OUTPUT COEFFICIENTS YN 
FOR A GIVEN SET OF INPUT COEFFICIENTS LFROM MO IN) 

COMPL EX DM'! 1P.1L YN' 1 1 i, h(1 1.5, 5), B>. 1 1, 5, 5) 

COrtriuri CuV HBnVTSC 5- 5 a- tBKP l 2* 5, 7>}> UuhP l l 25< 2> s, 5), PMGAIN<25. 
12. 5. 5). NLSYMC5. 5). NH(5, 5), RB-'tO. 5, 5), M0PPERC5. 5), N0PDER<5, 5), A, B 
COMMON HFLAG 

COMMON - GrtrP- WO .'POMN-' UMAX 

NLC'C FINDS THE OUTPUTS OF THE NONLINEAR ELEMENT 
IF LNFLAG.GE.3) MRITE <6, 30) . , ; „ ,, . 

ChuL NlOl <UH, ESKPP 1. i, N, NPP), UBKPT( 1, 1. H, NPP), PWGAINlI, i, N, NPP), 
INR'f'rSCU- NPP). NH^M- NPP), ABC1* N. NPP), NLSYM<N- NPP)> 

IF 1 NFLriG, GE, 3 > MRITE '6- 40) 
s iNLD=0. 0 
NK-NLSYI'UN. NPPV: 1 
HO 10 J-NK. UMAX 

GJW ChLCULhT £S THE VhLUE OF THE LINEAR ELEMENT FOR THE 
GIVEN VhLUE OF W 

_ ChLL Ij.’I'I CG. J- MO, H- l, N. NPP’G ES<i*H. NPP). MGRDERCN- NPPMIGRBERlN, No 
11 ’iF -.NFLhG.GE.3) WRITE <6, SCO 

1HE OUTPUT YN VaLL'E3_hPE FOUND BY MULTIPLYING THE NONLINEAR ELE 

outputs by the tprisfek function value of she linear element. 

VNC.D'G^DMCJ) 

I h (j r GT.b 1 GO TO 50 

10 continue 

IF > HFLhCL GE. 2" 1 MRITE (6-60) lDMvK), K=l, JMAN) 

IF lNFLaG.GE.3) WRITE Lb, 70; 

RETURN 
20 CONTINUE 

IF LNFLAG.&E.3) WRITE Lb, STj> 

STOP 


30 FORM!-: i ('20.1. 'ENTER GU I ’ ) 

40 FOPMhT (20X. 'PECNTEP OUT ’ ' 

50 FOPnfiT 20.'.. ’REENTER OUT' 

GO i Oft ,,1 iTT -.sh. 'Dr-i - .■ ilG’ 3 .S. .*. ■) 

7 0 rnpfiAT • TO' 1 . I F. -iyT- hoi ) 

SO FOPMnT <• -- '. 2QX. * NUMFJF 17, Oi- HhPNONIC TERMS EXCEEDS 5') 


C INVD'-'T 1 1 iE X MATRIX 
C 

SUBROUTINE NhTIIJV C..N. Y, IfcF NMhX) 

TilfiEMSIUN X' NMhX* HMA, - '..' ■ V' KMHr., k'lAH). jC'.25>. V’-‘2 ' 
COMMON I IFLhG 

IF * NFLmG. GE. 3 1 WRITE <6. rfl) 

Jt'p-u 

Lin 10 I =1. NMn:: 

DO 10 1=1. UMAX 
VI. J-XU. J. 

10 CONTINUE 

v.'n=i. 

CALL G jR f Y, NrirtX. Hfih'-:. H> ! f • *\? 0 . JC . v‘* 

GO TO 30 
20 MRITE (6,40) 

MR 1 1 L >'6, 50 * JC< 1 ) • V f ‘ i i 
1ER=2 

30 IF_' NFLnG.CE. 3’ WRITE *6,60.' 

PE I URN 
C 

40 FORMAT k 1 OX, ' OVERFLOW ' ) 

50 hOP.MnT 1 2 OX. £• G13. 7, i 01 1) '• 

60 FORMA T L40N. 'LEAVE MhTINV'* 

70 FOPMhT 1 4 OX, ' ENTER MhTINV') 

C 

END 



?nr in onooo 
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SUBROUTINE NLQC <DM, EbKPT- UBKPT, PWGAIN, NBMPTS- NH- A/ NSYM5 
DOUBLE; PRECISION Hi 

COMPLEX CH<11 >, ARGUE. APG- DMC11'». JMC115, SUM1, SUM2, SUM3 

DIMENSION EINC20015, JINDXB^QQ), KINDXTC5005, ISU<500> 
DIMENSIur! TB<S00>. XKt5005, bE1A<5005 

DIMENSION EBKPTC25* 25. UBKPT<25, 25, PUGAIN<25, 2), HflO) 
COMMON 'NL/' MPMX. I1EPNO, EPS /NOPT/ CH /QMUs NO /POMN/ JMhX 
COMMON Hf LAG 

APC,CNN. TZ5=CI-1PLX<0. 0, FLOAT'.NN5'M»IO*T2,> 

T>y J5=DELT+FLGAT<Jrl 5 1-0. 001 fDELT 
PI=3. 1 4 i 502653589793233 
PEP I OD=£ . t P T/ABS < MO > 

NSTP=HPH>:-1 

NpWG-NBKPTS+i 

DELT=PEPIOD/'NSTP 

IF >.NFLhG.G£. 35 WRITE <3. 44y> . 

1 1- • nfLhG.GE. 35 WRITE ■ G.'iBO.J NBKPTS 
DO 3 0 1-1,2 

IF kNFL AG.GE. i> WRITE 
IF NFLaG . GE • 3 j Nil I TE 
IP iNFLhG.GE.35 WRITE 
10 CONTINUE 

iF ^MFLA'j. GE. 55 NRIlE * S< 48 0-' 

IF <vNFLmG.GE.35 WRITE <6, 480' 

L-l 

NK 1 1 -■ = 0 - - 
BETA'- i I 1 = 0 . 


<3, 430 ) CEBkPKK, 15, K-b NBKPTS) 
f *• 480 > «: UBKPTCK, 1 5 , K--= 1 , NBKPTS > 
<3. 430; <PWGAIH I 'K, 15, K=l, UPWGJ 


'CN-:i5, I s 
HH, NSYM 


1.. jMaX'j 


J-l 

20 IF •’ J-NPfiN> 30,30.411 
30 T- TNlJ.‘ 

ElfK. I ■•=?'. T5 
J=.jh 
GO 10 20 
40 JFLiiG’J 
TBCij=0. 0 


Hv I ■ -- APRhY of hysteresis curve breakpoints 

Nil = NUMBER OF HYSTERESIS CURVE BREAKPOINTS 
-•SL~CUN i ROL FOR HVSTEPhSIS CURVE Stf tCTIuN 


SEircr INITThL NONLINETiRIIY curve 

JSW-l 

IF '.EIN f 2j.LT.EIH'*l '5 JSl.1^2 
IF Cin-I.LE. 0 5 JSW- 1 


ChLCI'LAIE ALL BFEhFPuINT TIMFS BY THE METHOD OF FALSE POSITION, 
WHICH IS ALSO \ MOWN AS THE SEChNT METHOD. 

DO 90 JJ=2- NBKPTS 
JJ2-JJ 

IF <*ElHCl5.L r.hBFPT' NBKpIS. JSW >5 GO TO 00 
JJ2=JJH 
JJi=JJ+l 
GO To 80 

50 CONTINUE : 

IF '1EINC1 i.GT. EBKPT <■ 1 ■ JSN55 GO TO SO 
JJi-JJ-i 
GO TO SO 

30 IF • T'EIN'15-EBTFTCJJ-I.JSl.O -.GT. 0. O '.AND. k<.ElN'.l5-EBKPT<JJ. JSW; 
15. LI. 0. ' • GO TO , 0 
GO TO 90 
TO JJ1=JJ 

80 XK''25-PMGAtri> JJi . JSW 5 

BETA- 15 -UUK.PT JJ2-1. JSW'-Fwn.-.UU.-JJl . .JSbKMEBkF’l‘<rjJ2-l . JSW . 

90 CONTINUE 

DO 250 nK-i.nSTP 
C DETERMINE EP&PKPOIN T TIMES 

K c W -0 
K-SGN-* !-l 

DU 250 JJ-1. NBKPTS 

IF <KSW.Fi'!. 1; GO TO 250 

If > ■ EIHu K5--FBKP i C J J. JSW5 « i KK+i -■ -FF’i- |-*T ■' JJ. JSW’ ’•> 1 00, 1 00 

1.250 
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no IF UJ.E0.HBKPT3) GO TO 120 

, , . - 1 *iV^ EI M c K " K ' " EBt: PT f ‘ J -H 1 • JSH*‘«=»<EIfKM*+i>EBKPT<JJ+l> » llu 

1 > lit 1 ’ 12 0 
110 kSGH=-l 

IF CitIrKKK>-EINChK+r*>.LT.O.> GO TO 250 
120 IF (ElfMTO.EGkEBiPTUJ, JSW.O GO TO 250 

K.SU-1 
II=II+i 

jinjd: :B<n>=jj*rsGH 
kiND>:Tdn=kr 

KSGH-+1 

ITER-i 

L=l+1 

:':i=T:kk'k> 

[ !2~1 ! J £ kK+ 1 ^ 

F i=E 111 CkK) -FBKPT < JJ. JSW> 

F2=EINCkK+l ,'-EBFPT(JJ. JSM ' 

IF CHFLhG.GE.3') WRITE (6, 430> 

IF (NFlhG.GE.3j WRITE (6.4r0,> k'k, JJ. II, XL N2* FI, F2, EIN(kK), EIN 
1 UCK+1 J ■ FBk'PT- JJ. JSl.D 

IF ivNFLHG.GE.iy WRITE 480> JSM 
IF 1 Fi.LU.h2,' STOP 
130 N3=::2-o:2-j:i>tF2-' r-2~Fi> 

IF (r2i.ECL0.02 GO TO 190 
F3=E ( '-..j > -EBr PT( J J,- JSt-i '■> 

IF I'ABSOil-Xll-EPS- 190.140.140 
140 IF < hBS(\3-N2)-EPS‘ ) 190. 150, 150 

it)0 It itl-tF3' i70- 160, ISO 

160 IF ■ hBS 1 FI ) . LT. EPS J GO 10 190 

170 F2=F3 

ITEP-HEP+i 

IF r iTtR.EO. ITERHO' GO TO 190 
GO TO 130 
ISO Fi=Fi 

xi=; 3 

iTER-ITEP+l 

IF ■ (TEF’.tii. ITEF’hO’ 1 GO TO 190 
GO TO 130 
190 TB'vL '=N3 

HE I’l.PMTNE WHICH CURVE 10 USE WHEN HVS1CKESIS UkEiskHUlHIS 
AFF. CROSSED 


i P • !'k,EO. i 1 GO TO 24 U 
IF (f)H.LF, 0 < GO TO 2-1 U 
IF 1 EIHkKK+i >.GT.H , 'iJ v GO IU 200 

351,1=1 

GO TO 240 

200 IF CtlFKH k+i ' . L i . rt* rirf • 1 GO TO 2i 0 
JEW-2 
GO TO 240 
21 0 Cum i Ihi it 

Nl=iHH-l> 

HO 220 1-1, HI 

220 IF k (EIN’J'k+i 1 . cE. h(I '■ 1 . hmm. f E 1 . L t . n f t +1 > > 

STOP 

230 IF •.EltW'K'.LE.AtlH -W--1 

IF ■FIrKi Y > . GE . A * 1 + 1 ' JSi.'-2 
240 CONTINUE 

IF dl-NSTP) 250.250.200 
250 CONTINUE 

ISiOF-ll-1 

260 TB d 1 j -PE.P I f iD+ 0 . OOHDELT 
C ORDER PLACES THE BREAKPOINTS IN THE CORRECT T I HE y-PUENCE 
CALL OPOEk (TB. JINB.-.B. RlNDXT- H • 


GO 


TO 230 


THIS DO LOOP KEEPS TRhCR OF WHICH PIECEWISE-L INEYiR GAIN 
GOES Li I iH EhlH BpEhePOImT llnE. 
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C 


C 


no 3io m::=2, ii 
hlF=4i: 1 

IF 1 1 > 270.300,300 

270 JM=IAB3< JINDXB<MX5 ) 

KX-KINDXT* MX' 1 

C . DETERMINE PROPER NONLINEARITY CURVE 
J5U=ISbkM>::' 

C MATCH GAIN HhD BREhKPOINT 

INC=0 

edot^eiikic :+i j~eim<k:o 

IF t < J3W.EQ. 2 ). hHI'. kJINBNEKliXJ.LT. 0)> IHC=1 

IF << JCtl.EO. 1 T.hND.CJINDXB' MX>.LT. OD ING=-1 

IF OEIiOT.LT. 0. >.HND. aUNDXBkilXJ.LT. 0. OANB. CNH.E9. 0>> INOl 

ii- < f edo r. gt.'O. :• . hub. •: jinu: :bcmx> . lt. o. t . and. <. nh.es. o > > inc=-i 

JX=JN+INC 

IF TEBOTT 290,420,200 
2S0 JX-JX-i 

290 BETA k ME .• =UBK PT C ..K J . • JU^-EET PT' JN, JSMUPWGAINC J.X+1, -JSWJ 

XK>' MX>”PI«IGAIHkJX+l, JSW> 

CO 10 3i0 

300 BET h ( MS 1 -BETh 1 1 ) 

XI kMX>=XK‘ 2.' 

310 CONTINUE 

IF vHFLAC.GE.3' 1 WRITE CB.4F0'- CEETACD, XK<IJ, 1=1, II) 

C 

C USt DFEAFPOINT TIMES hNB P’IECRWTSE-LINEaR GAlN INFORMATION TO ‘ 
C INTECRhTE TO GET hLL DM CD. 


no 410 p‘ri- 1 - jmax 

M=nn-i 

IH'MMV-=CMFLX< 0 . 0 . 0 . 0 - 
IF ai.EO .05 GO TO 330 
BO 320 1 = 3 . I STOP 

320 IMaitl ■ = 1 P; ^ MM BETA k I J i-kCENPCi iRGk- M. TBkl+1 D^-CENPCARG* -M. TB<I) 

1 > 1 • 

ipi- pui i M k Plri > +C MPlX k 0 . 0 , l .- k 2 . -rpl t F lOhT kri > ; > 

GO TO ^50 

330 DO 340 I=l> ISTOP 

^4 0 I pi v piri 1 = I n k mpF'+BETh i D * ■: TO • i + i 1 -TB * D j 

IMkMf-r-=IMCMM ■ 'PERIOD 
350 JM 1 Mil <-CMPL:X0 . U.uj. 0.< 

BO 4 00 I=i. I ST up 
SUMl=CMPLX<Ti. 0 . 0 . 0 ■ 

SUi 12 =CMPLXkO.O. 0 . 0 < 

SUM^- CMPLXk 0.0.0. Co 
SUMl“CN<l*<+kTBkI+D-TB.'DO 
DO 380 N=2, JMAX 

IF krt.tQ.M-i 1 GO ! 0 ?Af> 

1 > ' CHPLNk O^fwfen OAT oSd ‘ * M ‘ h TE * I+1> ’ ^™‘<RG<H M- 1, T Bk I > 
GO lO 370 

360 hPGUE=CN' N.' i <TB<I+l'.-TB f D^ 

370 SUM 2 =Si 'M 2 -.-hP.GUE 

jBU CONTINUE 

DO 390 Ib2. JMAX 

u . h , . , f^Mlf SUM3+C0N JGkCNkN'. 5 f- kCEi IPCARGk-N-M :•] TBO rl , ■ , -OXPkAPGk- 
_ . lN-N+1 . ! i fcf 1 1 -■ > "i fiPLFF 0. 0, WOfFLOHTk-N-MH 11 
4uu -'M 1 MMJ= JI-1' MMJ+k'SUMlTSI IM2+SI W +XLk I -Pi \ PERIOD 
BMkMM-fIMii-S-15+JMkl.iM. ' 

41 U CUM II HUE 

IF IhFUC 1 : 3 ? - ilSf ;’ c »' > 1 L 1 -- 16 -X 0 O WO ... 0.0 

RE l URN 
420 STOP 


430 FORMAT 1 T7, ’ Li. ' , Tia. ■* JJ' , T30- ' 1 1 ' . T42, ‘ l 'l * , r^i. ■ ■ '■’> t.'J ■pi , v fr - > 
^^. TQO. ’ EIN> NO . TS2, ’ EINO K-i D T1 04. ' EBKPf 1 jj'v'i “ 

■Wu FuRfiAT k2uX. ENTER NLOC* .■ 

1 2 OX, LEAVE Mi_OC‘ '■ 

55r! £S£' 1h I V"'?:7 BE Hi 71 ■■ 1F12.7, 10X. ’XL' = IF IP. 7. HTP 1 
480 FOPmhT ' 


END 
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c 


c 


£ 


c 


c 


SUBROUTINE FUNY FY, JLItf. YN. JFlhG, CM, UMAX, COST, fit) 

COMPLEX CHOU, Ylkll) 

COMMON HFLHG 

DIMENSION FYC25) • JFLRG<25j 
IF •.MFLAG.GE.3 - ' WRITE <£, 130-> 

DO 10 I-l.JMAX 

IF ( 'NFLhC.GT.3> WRITE <6. 90) I. YN<I V I, CN'O 
10 CONTINUE 

FV IS THE DIFFERENCE BETWEEN THE OUTPUT AND THE ASSUMED INPUT 

fs’< i >= i ;pehL':yN‘' i '> i-peaL'chcd v *■ 

DO 20 K-2. JLIM 
L=v2»kJ-2 

FY •: L ^ REhL C Yl KK) ) -PEALC CNC I : ) ) ) 

F'i v M > = *'hINrQ< YN •' F’ 1 )-rtI MAG <C NcK) ) 1 
20 CuH i I Nut 

COSTCO. 0 
DO 30 1-1, -JLJTi 

i0 CnST«Cus HUBS' FY<r- 1 

IF ■'PIF.El'.O'' RETURN 
INDEN=C. 

BO 4U LF-i, JLlM 

IF «. .IFLhGU:F) . EO. 0) IHDEX=INDEX+1 
FYOF- -FV'';'F+ INDEX) 

40 Cun 1 1 hi it 

L£it=JlIm- index 
KFLaG=0 

COMPUTE THE NORMALIZING FACTOR FOP THE OUTPU'I 
YNGR=0. 0 
BO 50 r-l/JMOX 

*>MOP=VH0F: j -^REhLuYNCK 1 1 j : *2 
S tiOF:=YlIOP+'.iiIMHiT-‘) i F.L j i t? 

50 CONTINUE 

YNUP-SORT 1 YNOP'> 

IF C 1, flt-OSj l-fl I i~l f'H 

CHECK FOFs ERROR WITHiN ShI I Si-', iL TORY LIMITS 

DO SU H-l.LIM 

It ■ i 'YiGS''h V> Yy '< > ’i ' ft iuf' J . G l . 0, 001 ' rtLAG-i 
EO CONTINUE 

IF 1 KFLwG. EC. 0) GO 10 70 
Rt l U C 'N 

70 WRITE C3. 100' 

WRITE 13fN> 

RETURN 1 

SO WRITfc C6. 1 1 0.< 

STOP 


90 FOPMhT czOX. 
i.' 

100 FOPMhT '.a 
110 FORMAT -- 20 X. 
120 FORMAT -.2 OX. 
130 pOpmaT '.2 Ox. 


•YiU', 12. 0=',2'F12. 7, 5X),?::. 

SOX, ' ERROR LIMIT SATISFIED ' v 
hLL OUTPUTS ARE ZERO' 

•' ENTER FUNY' *> 

AC 1 tP I KrEii.E Vs-.LOt'S hPk ' ’> 


' L‘N ’ M2. 




la./.sx* 


ENP 



non on n no n non o o o non non o o, 
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SUBROUTINE READ COLLECTS INPUT DATA. 

INPUT BfiTA CONSISTS OF INITIAL GUESSES AT INPUT <CN> 
COEFFICIENTS AND FREQUENCY (WO) FOR HARMONIC BALANCE, 
BREAKPOINTS AND SLOPES WHICH DEFINE THE NONLINEARITY, 
AND COEFFICIENTS OF THE LINEAR ELEMENT POLYNOMIALS. 
ALSO INCLUDED ARE THE NUMBER OF PARALLEL PATHS AND THE 
NUMBER aND FORM OF NONLINEAR ELEMENTS IN EACH PATH. • 


NOHAR » NUMBER OF HARMONICS < EXCLUDING DC ) 

JhhX - NOHAR + i.= NUMBER UF HARMONICS INCLUDING DC 
DELIA IS' THE PERTURBATION FOR FINDING J 
JD IS THE CURVE SELECTOR 

AB = COMMON BREAKPOINTS UN HYSTERESIS CURVE 
NLSYM=0 FOR ASYM ML. 1 FOR SYM NL 
ITERNO= HAN. NO. OF ITERATIONS FOR NL OUTPUT 
£F'S= SECANT METHOD STOPPING LIMIT 

hPMP.= NO. OF TIME POINTS USED TO DESCRIBE INPUT FUNCTION 
FF= MINIMUM STEP SIZE FOP TERMINATION 
NORDtP. M0PDtR= ORDER OF DENOM. AND NUMER; OF LIN. EL. 
SIZE DC f ERMINES THE LARGEST CHANGE IN CN WHICH ChN OCCUR. 


SUBROUTINE READ 

COMPLEX CHI CM’ x Ai.ll, 5- 5), B(li.5,5> 

DintNSION HAP.C22). hACIIx PB'ii* 

•COMMON NFL AG 

common . Maine. - cni, delta. ee,size. hohar 

L'UMM'JN 'OmF' iilu -puMN/ UMAX -ML/ NPNX» J.TEPNQ, ti-S 
COMMON -'OPT/ NPj-NNL 1 ' 5-’ 

COMMON /QT/ NBRPTSC5. SX.EBKPT'ZG, 2, 5- 5), UBKPTC25, 2, 5, 5), PUGAINC25, 
.12, 5. 3), HLSYfl’XL OX. NH(5, 5 *. HBklO. 5. 5), H0RPER(5, 5x H0F’DFR(5, 5X A, B 
WRITE u5, ll.O’ 1 

READ (5, ISO) DELTA, NPMX, ITEPNO, EPS. EE, SIZE, NFLAG 
READ o, 140X HOHAR, HP 
IF f DELTA. EO. 0. .< DELTh=1 . E-OS 
.IF <NFMN.EQ.O> NPMN-SOl 
IF (It EPNU. bU. 0) I (EPnO=i 00 
IF * EPS. EO. 0. ) EPS=1.E“G5 
IF kSIZE. EO. O'. > SIZE=-0. 1 
"IF (NP.tU.0'> HP=l 
IF ■'EE.EU.O) EE=l.E-08 
IF CNFLAG.GE.3’ WRITE i.F. 120-> 

WRITE i6, 1F0) 

WRITE CF.180'* DELTA, NPMX, 1TERNO, EPS 
WRITE (S. 400; 

WRITE > 6. 3?G * HP LAG. E£. SIZE 

JMAX=NOHhR+1 

IriD=2*N0HAP 

i.jR 1 1 e <e. 150- JriA>: 

READ <5. 140; WO 

READ (5,140) (HaR(Jx J-l. IND) 

READ (.5, 14 O’ 1 CHNLc lx 1=1, HP) 

WRITE /». 170) NOHAR 
WRITE (b, 130 J WO 

WRi ! t (6, £0U) '.NP. i.HNLCI.'. 1=1- HP. 1 ) 

WRITE (b, 210' 1 
WRITE ib. 230' HhR(1) 

WPITE (b- 210 ' 

-WRITE • b. 240' HaRi>2) 

IF (NOHAR. EO. 1) GO TO 20 
DO 1 0 i i =2 , NQHaR 

10 WRITE iS, 2501 IT, HARi24lT-i:'.HAR(2tir) 

20 DO 80 KP=1 • HP 
Hl=Hnl- kp"' 



92 


DO 80 K=l, NL 

READ <5, 140) NBKPTSCK, KP), NLSYM<K,'KP> 

• WRITE Cb, £60) K, kP, NBkPTSCK, KP) 

IF <NLSYMkK,k'P).EC!.0) WRITE (6, 270) K, KP 
IF <NLSYM<K, KP) . EGL 1') WRITE <6, 230) K, KP 
IC'P=NBKPTS(R, RPV+1 
Kk=NBKPTS<K, KP) 

30 READ <5, 140) JD 

IF CJD.EO.'O) GO TO 40 
WRI lb <6, 290) JD, K, KP 
READ r5, 140) fEBKPTCI, JD, K,KP), 1=1, KK> 

READ *'5. 140) <UPKPT<I, JD, K* KP), 1=1, KK) 

PEhU '5, 140 - <PWGAIHa, JD, K, KP), 1=1, IOP) 

WRITE <6,300) (EBKPTkl. JD, K, KP)/ 1=1, KK) 

WRITE -t, 31 0) (UBKPTa, JD, K, KP), I5I, KK) _ 

WRITE c 6, 320) kPWGAlf'R I, JD, K, KP), 1=1, IRP)' 

GO TO 30 

40 READ <5- 140) NH<K. KP) 

NK’=HH<R. KP) 

IF (NHCDKP'-.EO, 0"> GO TO 50 
WRITE kb, 330) NH(K.KP) 

Read (5, i4D) chBU. k. KP), :I=i, NX) - 
WRITE -.6,340) <AB<I» K* KP' 1 , 1=1- NX) 

50 CONTINUE 

PEA 0 '5, 140- MOFBER- K- KP), NOPBER<K» KP> 

WRITE <6, 350) MORDEPCK. kP), NGRBER<K,.KP.'-. K. KP 
MTERMS=MOPBERkK. KP>+1 . . ‘ 

HTERnS=NOADtPkK, kP)+1 

READ <5, 140) <hA- X>, J~1 . MTERMS’ 1 

■F-EhD <5,140) <DB<J). J=l, NTERMS) 

WRI 1 1 <6 , 360) Chmi J - . J= i, M I ERMS) 

WRITE kb. 370) <BBkJ,, J=l, NTEPMS) 

DO 60 J=l. MTEPMS 

60 A<J. K- kP ) =CMPLX< Ah < J-> ,‘ u . 0) 

DO 70 J=l. NTERMS' ’ 

70 B< J. K- kP)=CNFL::>BBO), u. O'- 

SO CON 1 1 Mi iE 

CNI<r)=CMPLX<HAPki), 0. O'- 
CNIi2 -CHPLN-.O. 0, <HhR<2) V-2. ) 

IF 'NOHaR. E l 1 , i) GO IG too 
DO 90 H=3. JMaX 
M=24N 

90 CHI <N)=Cf-rLK< < mP- h--*.' '-'S. • -riMR'. M-2 V-:-?. ) ) ) 

100 CONTINUE 

URITE_<b-410) 

If- -. NhLAG.Gt.'i) MPITt <6- 130) 

RETURN 

110 FORMAT <lHl,,v'/''i0i:, !£'<' + JJiHPUT DATA MCK' +')) 

120 FORMAT <2 OX- ’ ENTER READ - ’ '• 

130 FORMAT <20::. -’LEhYE READ") 

DO FORMAT O 

150 FORMAT <, 10X. 'JNrtX = ' ■ 12) 

160 FORMAT <///.•/, 5N, 'DELTA', T15, 'HO. POINTS', T29, 'NL ITER', T48, 'SEC. 
1STOP LIM' > 

170 FOP.MhT -V^ICX, 'NO. OF HARMONICS =', ID 
180 FORMAT <E10.5, £110. 3E10.5. 15)’ 

1 9 0 1- QRMi A •* ' - 1 OX ■ I N I T:I AL U =', Gl 7) 

£00 FuRMsiT <."J OX. ' HUMBER OF PaIHS = \ I U, 3X, ' NUMBER OF',' NONUNEARIT 
1TES FOR PhTH = ',5I4- 
210 FOkMhT >V.'2K{. ' INITIAL INPUT VALUES' ) 

220 FOF’MuT * -VT25. 'HhRMOHIC VhLUES', 745, 'MAR. NO. ', T57- 'COSINES', T72, ' 
1SINESD 

230 FOkmAi '.-•‘''20J. 'EC =‘,Gl£.7'' 

240 FOPIIhT ■ '2 ON, 'FUNDAMENTAL COMF, ', T70, G12.7) 

250 FORMaT k/T47, l£. T55, G12. 7, T70. GJ 2. 7) 

260 FORMAT •‘.'✓I X. 'NO. OF BREhK POINTS IN NONLINEARITY', 12. 1 OF PATH', 12 
1. ' IS '.ID 

270 FORMAT C.-.'l.X. 'NONLINEARITY’, 12. ' OF PATH'. 12, ' IS ASSYNETPICAL ' > 
280 FOPMmT '.ssW, 'NONLINEARITY' - I2< ' OF PhTH'.JB..' IS SYMhTRICAL') 


CL 
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290 FORMAT (MX, 'BREAKPOINT VALUES FOR CURVE M2, ' OF NONLINEARITYM 
12- ' IN PhTH' ■ 12) 

3U0 FOPMhT (',5X, 'HORIZONTAL COORDINATES', T30, S(E1 0.5, 5X)> 

310 FORMAT C'.'SX, 'VERTICAL COORDINATES'. T3G, 8<E10.5, 5X>> 

320 FORMAT t"3X, 'SLOPE OF SEGMENT TO LEFT', T30, 8(E1 0.5, 5X)> 

330 FOPMhT '1'fiX. 'NUMBER OF COMMON POINTS ON H7S. CURVES -M2) 
i4U FORMA l 'W2X, 'HORIZONTAL COORDINATES OF COMMON POINTS' , 3N- 8(E1 0. 5, 
14X> ) • - 

350 FOPMhT f/"2X, 12, ' ORDER NUMERATOR, M2, ' ORDER DENOM.',' IN NONLIN 
iErtPiTY', 12, ' IN PATH', 12) 

300 FORMAT t^SX, 'COEF.OF NUM. TERMS (HIGHEST FIRST)', 8<5X, E10.5)) 

370 FORMAT (MX. 'CGEF.OF DEN. TERMS (HIGHEST FIRST V, 3<5X. El 0.5)) 

380 FORMAT 'lX. E10.5, 2110, T50, E10.5) 

390 FORMAT (IX) 110, T25, El 0.5. T50, El 0.5) 

400 FORMAT K' '5X. 'PRINT FLAG' . T25, 'MIN. STEP SIZE', T45, ' INIT. STEP') 


C 

c 

c 


c 

c 


Finn THE TPhNSFEP magnitude and phase of the lineup system describe 

BY A hNB B AT THE FREQUENCY l.*0 


SUBROUTINE GJM (G. J, WO. A, B- MOPDER, NOPDER) 
COMPLEX G. Z. P, H( 1 1 ), Oa 1 ) 


COMMON r If- LAG 
IF ir!FLAG.GE.3T 


WRITE (8.50 % ' 


Z=hU> 

IF (MOPDEP.EQ.O) GO TO 20 

nn i n M?=1 . HfiKTlFR 

1 0 2-Z+CMPLXC 0. 0. FLOAT >' JL> *U0 '+h(MZ+1 .* 

20 CONTINUE 
P=B'.l.' 

IF • NOPDEP.EO. 0.' GO lO 4 0 
DO 30 NP=1, NOPDER 

30 P=PtCMPL/XO. 0, FLO aT< JL j * WO ) •! B <. NP+ 1 T 

40 CONI I Nut 

IF ^CABS(P).EO.O.O) p=p+CMPLX‘. 1 . E-30. 0. ) 
G=Z P 

if ■.peai/g).f:o,o.o> g -G- f c mpl: ; c i . e- 3 o * o . ) 

IF (NFLhP.GE. 2T WRITE ^6, 60> G 
IF f NFLrtG.GE.3> WRITE • 6. 70'- 
RETURN 


50 FORMAT ■ SOX. ' ENTER GJW”> 

60 FORMhT -mX, 'GU3EAL,' -'.G15.7, 10X, 'gumag:* 
70 FOFMAT *.2 OX. 'lEAVE G J.F .• 


- G15. 7) 


C 


END 



GO 
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c 


c 


c 


c 


c 


DETERMINE THE CHhNGE TO BE MaDE in the assumed inputs 

SUBROUTINE DELC kDELCN, flJIHV, JLIM, FY, KBIG, JFLAG, DUMMY, SI2E) 
DIMENSION DElCN< 25), fiJINV'25. 25>. FYk25>, JFLAGk25> 

COMMON NFLAG 

IF CNFLhG.GE.3) WRITE <6, S0> 

DO 10 1=1, KBIG 
ID CONTINUE 
DO 20 K=l, JLJM 
DELCN<K>=0. 0 
DO 20 J=i, JLIM 

BELCH k K ^ < A J 1 N V k K • *FYU»+DELCN(K> 

20 CONTINUE 

ADJUST BELCH TO ACCOUNT R3R ZERO ROWS IN J MATRIX 
NISH I FT- JL I M “KB I G 
IF CNSHIi-T.tO. 0> GO 1 0 -10 
DO 30 JVrtR=JLIM. l.-l 

IF k JFLNG k J V AP } , ECU 1 T BELCH C J VAR > =DELCH C J V AP-HSH I FT *» 

Ih t JFLflG' JY hRj.EG'O) BELCH < J V AR >= 0. 0 
IF UFLnGUVAR>.EO.OJ NBHIFT=-NSHIFT-t 
30 CONTINUE 
40 CuMl tHUFI 

DETERMINE THE LhRGEST CHANGE 
BUMMY^O. 0 
DO JO \ i=i, JLIM 

IF CaBSC BELCH* ^ J, GT. DUMMY J BUi \\ IY~i iBS ^ DELC N k M > } 

50 CONTINUE 

wpijL c 6 « i o o 1 dummy 

IF THE LARGEST CHANGE IS GREATER THAN "SI2E‘, SChLE DOWN 
IF kDUMMY.LT.SIZE> GO TO 70 
DO 60 KZ"1 j AIM 

DELC N k K2> -DELCN k K2 > **S 1 Z£s DUMMY 
60 CONTINUE 
BUMilV-SlEE 
MRI1F kb. 110 « DUMMY 
70 CONTINUE 

IF ’Ni-lhG. GE. 3 ' WRITE k6»S0.« 

RETURN 


8 0 FORM A t 1 2GH- EN ! t P DELC * > 

SO FORMA I 'gQK» ' LEAVE BElC* ? 

format <. 20 ::, ••max. delta comploed =%gi 2 .?> 

llU FOPMiiT >.1H+, TbOj •’Mn::. DELTA ALLOWED ='/ G12.7; 


END 
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APPENDIX C 


SAMPLE OUTPUT OF HBA 
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***** 4 ***«JUPUT 


data *+♦**♦*+* 


• XoBcbias rtC> - W TS » 0 XTcK 

PKIHT FLAb Ml*-!. i>Trp S12S. 

0 •i0\j^0-0T 

UMAX = •* 

KO0 bK h ^ mO . nIC ^ -3 
INITIAL * = 2.50y000 


SEC 


aSJfl-W' 


INiT. i>Tc° 

.lyOQO+OO 


NU w c=.H Or PATHS =1 

>|US;»EK OF i'i^/%Ll ,1 *EA^iTlt.5 FOR PaTh - 3 

initial input values „ 
oc = *000^000 

HAkr.-0.NlC ^mLU^S HArt • NO* COSINES 

fundamental comp, 

2 .0000000. 

3 .0000000 

NO « 0- aRsAKPOXN T^> IN NO^LlNuARj T f 1 OF PATH X IS 3 
\u%LI.£aRITY 1 OF PATn 1 IS »YWteTRlCAL 

b*t*AP0INT VAuJ£S FOR CUK<E 1 OF .<OU-INEaRXTY 1 IN PATH 1 

Hv.RIZO'11 AL COUPpJ.NATcS “*10000+06 -.20000 + 0l .20000+01 

vertical coordinates -.30^00+01 — 3 ogoo + qi .30000+01 

SuOHE OF SEVENS TO LE^T .GO&vO .00000 .15000+01 

I (rOER f,UiV,L c <AT^?tf 2 G*S£R OENOX. IN NCNLINEaRIT* 1 IN PATH l 

» i » 

C'tF * CF nUm* TeRinSMISH-ST FIRST) .50000+01 .20000+02 

CoEF.OF 3EN. tEK*!>(HX3r.&bT MRSD .10000+01 .30000+01 

NO. 0- |5A£AKPU1||T» IN NOolI^ARITT 2 OF PATH 1 IS * 


.jw Ll JEAHlTY * OF PAT* 1 iS i>t*' -iRi^Al 

£>Rcan^OI‘mT \ZAj-ueS Fo*R CuHvE 1 OF NGNL1.N£AkITY 2 In PATH 1 


ORIGINAL PAGE IS 
OF POOR'QyALITY 


SINES 

1,000000 

.ooooooo 

tOoOOODO 


• OOoOo 


.50000+01 



H Ja1£v‘‘JTAl C&^CJlNAT-S “♦lo < iOQ+uE> -*20000+00 

VERTICAL uOO^^lNATES -.1500'J+Jl -*1501)3 + 01 

SL'J^tL OF TO L E^T *OtWO *OOOCG 

ja-ACPOIM V*^£S rC-H CJ*%'E * OF ^LINEARITY 2 IJ PATH 1 
M^UC’lThL cOOH;I JATS.S -.lOCJG+06 -*200°0+u0 

vupTi cal coordinates -.ibouj+^i — .isooc+oi 

5L0-L Or Sw'uM TO L£r f -0C--0 •COOOfr 


•acooo+cc 

— IbOOO + Ol 
•00003 

-*20000+00 
•15000+01 
• lOCO 0+21 


TJvi^ Or LO.,vo.\ POINTS ON nYS.CuT/£S = 2 
M-RI^C.'JAL C^^rC‘ 1..AT^S Or Civics PORTS' -.20000 + 00 .20000+00 

2 0~<Ca* 3 O'jE* b*Xl'U AN NONLl.NiARUt 2 I?4 PATH 1 

Cc£r*0r NU* 1 - f '.b IH;oh~ 5T *10QCC+0l .30000 + 01 

CtiF.Or 3: u T t K ^(HiSH-ST rI*Rsn .IQCOO+Ol ,70000+01 

•^ f C- or'EAn^AjJ^ ri NO*LIn~AHiTY 3 Qf- PaTH l IS <5 
.<* uI^A'XTT i o* PAT.-, 1 IS sY lJ *. 5.TRJCAL 


•2UqOO + 00 

•l^OOO+Ol 

•10000+21 

• 20000+00 
• 1^000+0 i 
* OOOOQ 


.50000+01 

*15000+02 


vIjm^C IM V*u.£b CV<*5 1 OF lONUlNEARiTY 3 1*. PATH 1 

Cv^l.iAT-S -.lC‘)00+^6 ~*3OC0J+Ol *30000+01 

V'-.OICi'L ww^::^KAr:S '*25000+^1 -.25000+01 *25000+01 

vLO-^c. OF TO u£r T .CCCOO .'JCQOO .03333+00 .OOOOO 

; / i C'JEH IN MONU^-A^HT 3 IN path 1 

C^Er.Of* N‘J *\ * ) £rt* o t HI SH—ST “IRSD *10000+02 

C-Lr.Cr 3E'<* .3< nirrftf.br J"IRS?) 


•OOO9O 


.COOoO 


,12000+02 


loooo+oi 


^OOOC+Ol 



****** 4 ****** + 


results 


*************** 


*♦* 


COST 1 

♦ ^t-LTA 


- *10 


MmX • ^4t 1AL P^hVjrtOATlON - *1000000-05 
wOOOO+U COST 2 = .4697641 _ *** 
z.u = . 829 oC 23 MAX • DELTA ALLOrt=.D = . 


xqoooqo+oo 


Nc.ft wO = 2 *oOOOOO 

Ns. a DC COMPONENT = *0000 000 

NLrt FUNDAMENTAL SINE COMPONENT r .9506131 

Nc,; CONSONANTS FOR HARMONIC NO* 2* COSINE = 1.999915 

N“ r . COMPONENTS FOR HARMONIC HO* 3* COSINE = 2.004194 


SINE r -*1725285*03 
SINE = —9212824-02 


*** 

■I A A * wiLt*' 


^ P«3^f?AL 8 *[fiVb3^TlJ5 = 


• 43^7041 
'VJTtLi = .3110745 


20 QU 000 

C0& VAX. "DELT A 2 ALlO^LO^^ .2000000 


NZ* *0 = 2 * bOO 000 

Nift DC COMPONENT = ♦ 0000000 

n£;\ Fa\DA.w£NiAL 51 NE COMPONENT - • .8606977 

Ndrt COMPONENTS p^R PARSONIC NO* 2* COSJNE r 1.9Q9630 

»N£ rt CO 4PONt‘J I S FOR HARMONIC NO* 3» C0SlN£ = 2.012131 


SINE r —1327442-02 

sine - -* 2120595-01 


*** CUaT 1 
’» A X * 


**; R«. = . 4 C 0 C-O 0 O 


— • 4 i* T S 4 o 2 

COMPUTE = . 645/211 


'COST 2 = . 313614 b *** 

VAX. DELTA ALLOWED S .4000000 


N-* r.O = 3-200000 


Nt/ ( DC COvWJNEmT = .0000000 

*l~„ FU !TA U SINE COMPONENT = .6725511 

N- CO.MPO.NcNtS FOR PAR IONIC NO* 2? COSINE = 2.000208 

,L ( , CON-ON'.vrS FOR ‘■'ARMONIC NO* 3? COSINE = 2*022547 


SINE = *1514935-02 

SINE = -*2949087-01 


00 



*** RUN NJ’-MER 4 F? LLOtfS **♦ 

V,AA. PARTIAL PERTURBATION = .8000000 

*** C^ST 1 = 0136140 COST 2 - .13*5140 *** 

.VjAa* u-lTA CO'VuiO = .3559791 

Nt* liO r j 053979 

Hdit DC COMPONENT = .0000000 

NS,v FUND A -TENJ SINE COMPONENT = *5129287 " 

N£* COMPONENTS FOR HARMONIC NO* 2* C05JNE = 2,000385 

NEa COMPONENTS FOR HARMONIC NO* 3* COSINE = 2, 018022 


*** RUN NU \6C* 5 FOLLOWS *t* 

t Nea MA*. PA < 1 IAL PE^TJRBATIQN = .71195*3 

*** COST 1 = •13*614* * u CO£t 2 r *962554>2-U2***^ 

m>U. uuLTa computed = , 1932573-01 

Ne* WO = 3*575306 

Me*; DC CO 'POTENT = *0000000 

NErt FUNDAMENTAL SINE COMPONENT = ,5019447 

Ne„ COMPONENTS FOR HARMONIC NO* 2* .COSime = 2.0Q0006 
CO •i^ONLfiTS FOR HARMONIC NO* 3* CoSInE s 2,016803 


Sine = — 5514378-03 
SINE = -*2236249-01 


SINE = -*9780168-05 
SINE s -* 2183689-01 


r** CUST 1 


*** RUN NUVOC.R 6 FOLLOWS *** 
l'C„ M A A , P/»R l IAL PERTURBATION = *3W65346"0l 

.9o2555d-02 CO^T 2 = .7602895-04*** 


s°.*V i-IMiT SATlShI£9 

aCCePUJLL values are 

OC = .oooucoo 

FU.OAMenTAu SINE CO’-^OvEnT z .5019^4 7 

COMPONENTS For HARMONIC NO* 21 COSINE = *1220431-04 SINE = -*9780168-05 

CJM^lNe'ITs -C* U ARV,C iic NO* 3i COSINE = *3360616-01 Si ME = -*2!038fP-Ol 


<o 

<o 


L i v I T CYCLe FPEO^ENCr = 3.875306 


RADIANS/SEC 
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